A simple two-dimensional parameterisation for Flux Footprint Prediction (FFP)

Abstract

Flux footprint models are often used for interpretation of flux tower measurements, to estimate position and size of surface source areas, and the relative contribution of passive scalar sources to measured fluxes. Accurate knowledge of footprints is of crucial importance for any upscaling exercises from single site flux measurements to local or regional scale. Hence, footprint models are ultimately also of considerable importance for improved greenhouse gas budgeting. With increasing numbers of flux towers within large monitoring networks such as FluxNet, ICOS (Integrated Carbon Observation System), NEON (National Ecological Observatory Network), or AmeriFlux, and with increasing temporal range of observations from such towers (of the order of decades) and availability of airborne flux measurements, there has been an increasing demand for reliable footprint estimation. Even though several sophisticated footprint models have been developed in recent years, most are still not suitable for application to long time series, due to their high computational demands. Existing fast footprint models, on the other hand, are based on surface layer theory and hence are of restricted validity for real-case applications.

To remedy such shortcomings, we present the two-dimensional parameterisation for Flux Footprint Prediction (FFP), based on a novel scaling approach for the crosswind distribution of the flux footprint and on an improved version of the footprint parameterisation of . Compared to the latter, FFP now provides not only the extent but also the width and shape of footprint estimates, and explicit consideration of the effects of the surface roughness length. The footprint parameterisation has been developed and evaluated using simulations of the backward Lagrangian stochastic particle dispersion model LPDM-B . Like LPDM-B, the parameterisation is valid for a broad range of boundary layer conditions and measurement heights over the entire planetary boundary layer. Thus, it can provide footprint estimates for a wide range of real-case applications.

The new footprint parameterisation requires input that can be easily determined from, for example, flux tower measurements or airborne flux data. FFP can be applied to data of long-term monitoring programmes as well as be used for quick footprint estimates in the field, or for designing new sites.

Introduction

Flux footprint models are used to describe the spatial extent and position of the surface area that is contributing to a turbulent flux measurement at a specific point in time, for specific atmospheric conditions and surface characteristics. They are hence very important tools when it comes to interpretation of flux measurements of passive scalars, such as the greenhouse gases carbon dioxide (${\mathrm{CO}}_{\mathrm{2}}$), water vapour (${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$), or methane (${\mathrm{CH}}_{\mathrm{4}}$).

In recent years, the application of footprint models has become a standard task in analysis of measurements from flux towers e.g. or airborne flux measurements of inert (i.e. long-lived) greenhouse gases. Information about the sink and source location is even more crucial for towers in heterogeneous or disturbed landscapes e.g., a topical study area in recent times. showed that the use of the footprint concept has been increasing exponentially since 1972. Since then, footprint models have been used for planning and design of new flux towers and have been applied to most flux tower observations around the globe, to support the interpretation of such measurements.

Long-term and short-term flux observations are exposed to widely varying atmospheric conditions and their interpretation therefore involves an enormous amount of footprint calculations. Despite the widespread use of footprint models, the selection of a suitable model still poses a major challenge. Complex footprint models based on large-eddy simulation LES; e.g., Eulerian models of higher-order turbulence closure e.g., Lagrangian stochastic particle dispersion LPD; e.g., or a combination of LES and LPD e.g. can, to a certain degree, offer the ability to resolve complex flow structures (e.g. flow over a forest edge or a street canyon) and surface heterogeneity. However, to date they are laborious to run, still highly CPU intensive, and hence can in practice be applied only for case studies over selected hours or days, and most are constrained to a narrow range of atmospheric conditions. These complex models are not suited for dealing with the vast increase of long-term flux tower data or the more and more frequent airborne flux measurements. Instead, quick footprint estimates are needed, models that can deal with large amounts of input data, for example several years of half-hourly data points at several observational levels at multiple locations. For these reasons, analytical footprint models are often used as a compromise e.g.. These models are simple and fast, but their validity is often constrained to ranges of receptor heights and boundary layer conditions that are much more restricted than those commonly observed.

Existing footprint modelling studies offer the potential for simple parameterisations as, for example, proposed by , , or . The primary drawback of these parameterisations is their limitation to a particular turbulence scaling domain (often surface layer scaling), or to a limited range of stratifications. Conversely, measurement programmes of even just a few days are regularly exposed to conditions spanning several turbulence scaling domains.

To fill this gap, introduced a footprint parameterisation based on a fit to scaled footprint estimates derived from the Lagrangian stochastic particle dispersion model, LPDM-B . LPDM-B is one of very few LPD footprint models valid for a wide range of boundary layer stratifications and receptor heights. Likewise, the parameterisation of LPDM-B is also valid for outside surface layer conditions and for non-Gaussian turbulence, as for example for the convective boundary layer. However, the parameterisation of comprises only the crosswind-integrated footprint, i.e. it describes the footprint function's upwind extent but not its width.

Recently, footprint model outputs have frequently been combined with surface information, such as remote sensing data e.g.. As remote sensing data are increasingly available in high spatial resolution, a footprint often covers more than 1 pixel of remote sensing data; hence, there is a need for information on the crosswind spread of the footprint. Similar to the flux-source area model, FSAM , the footprint model of includes dispersion in crosswind direction. provided a crosswind extension of the footprint model of . However, all these models are of limited validity restricted to measurements close to the surface. For more details on their validity and restrictions, and for a comprehensive review on existing footprint techniques and approaches, the reader is referred to , , or .

This study addresses the issues and shortcomings mentioned above. We present the new parameterisation for Flux Footprint Prediction (FFP), with improved footprint predictions for elevated measurement heights in stable stratifications. The influence of the surface roughness has been implemented into the scaling approach explicitly. Further and most importantly, the new parameterisation also describes the crosswind spread of the footprint and hence, it is suitable for many practical applications. Like all footprint models that do not simulate the full time- and space-explicit flow, FFP implicitly assumes stationarity over the eddy-covariance integration period (typically 30 min) and horizontal homogeneity of the flow (but not of the scalar source/sink distribution). As in , the new parameterisation is based on a scaling approach of flux footprint results of the thoroughly tested Lagrangian footprint model LPDM-B . Its most important scaling variables are readily available through common turbulence measurements that are typically performed at flux tower sites. The code of FFP can be obtained in several platform-independent programming languages at www.footprint.kljun.net.

Footprint data set

Mathematically, the flux footprint, $f$, is the transfer function between sources or sinks of passive scalars at the surface, $Q_c$, and the turbulent flux, $F_c$, measured at a receptor at height $z_{\mathrm{m}}$ e.g.. We define a local footprint coordinate system, where the receptor is mounted above the origin ($\mathrm{0},\mathrm{0}$) and positive $x$ indicates upwind distance, such that $$F_c(\mathrm{0},\mathrm{0},z_{\mathrm{m}})=\underset{\mathrm{\Re }}{\int }{Q}_c(x,y)f(x,y)\mathrm{d}x\mathrm{d}y,$$ where $\mathrm{\Re }$ denotes the integration domain. As the footprint function is always specific to a given measurement height, the vertical reference in $f$ is neglected, for simplicity. It follows that the footprint function is proportional to the flux increment arising from a single unit point source or sink, $Q_{\text{u}}$, i.e. $$f(x,y)={\displaystyle \frac{F_c(\mathrm{0},\mathrm{0},z_{\mathrm{m}})}{Q_{\text{u}}(x,y)}}.$$ As $F_c$ is a flux density (per unit area) and $Q_{\text{u}}$ is a source or sink integrated over a unit area, the two-dimensional footprint function has the dimension of (1/area). Assuming that crosswind turbulent dispersion can be treated independently from vertical or streamwise transport, the footprint function can be expressed in terms of a crosswind-integrated footprint, $\stackrel{\mathrm{‾}}{f^y}$, and a crosswind dispersion function, $D_y$, see, e.g. $$f(x,y)=\stackrel{\mathrm{‾}}{f^y}\left(x\right)D_y.$$ In the following sections, we present a scaling approach and a parameterisation for the derivation of $\stackrel{\mathrm{‾}}{f^y}$ and $D_y$, with the aim of simple and accessible estimation of $f(x,y)$.

Derivation and evaluation of the footprint parameterisation are based on footprint calculations using LPDM-B . LPDM-B is a footprint model of the Lagrangian stochastic particle dispersion type, with three-dimensional dispersion of inert particles as described by and . LPDM-B fulfils the well-mixed condition and is valid for stable, neutral, and convective boundary layer stratifications, assuming stationary flow conditions. It has been shown to reproduce wind tunnel simulations very well . LPDM-B reflects particles fully elastically at the surface and at the top of the planetary boundary layer and tracks particles backward in time, from the receptor location to the source/sink, at the surface i.e. particle touchdown location; seefor details. Hence, footprint calculation can be based on all computed particle tracks directly without the need for coordinate transformation. We refer to , , and for details in the formulation and evaluation of the model.

Compared to the original parameterisation of , we have increased the parameter space for LPDM-B simulations especially for stable boundary layer conditions. We also increased the covered range of roughness lengths, $z_{\mathrm{0}}$, to include roughness lengths that may be found over sparse forest canopies. For the parameterisation, a total of 200 simulations were run with LPDM-B for measurement heights between 1 and 1000 m and boundary layer conditions from strongly convective, neutral, to strongly stable. With that, the simulations span a range of stability regimes, namely, the surface layer, local scaling layer, z-less scaling layer, neutral layer, the free convection layer, and the mixed layer cf. . Table  gives an overview of the parameter space for the simulated scenarios. We use standard definitions for the friction velocity, $u_{*}$, the convective velocity scale, $w_{*}$, and the Obukhov length, $L$ see, e.g.and Appendix  for details and for the definition of $L$. Each scenario was run for the whole set of roughness lengths and for all listed measurement heights. Note that we define the measurement height as $z_{\mathrm{m}}=z_{\mathrm{receptor}}-z_{\mathrm{d}}$, where $z_{\mathrm{receptor}}$ is the height of the receptor above ground and $z_{\mathrm{d}}$ is the zero-plane displacement height. 47 cases with measurement height within the roughness sublayer ($z_{*}$) were excluded from later analysis ($z_{*}\approx n{h}_{\mathrm{rs}}$; where commonly $\mathrm{2}\le n\le \mathrm{5}$ (), and $h_{\mathrm{rs}}$ is the mean height of the roughness elements, approximated by $h_{\mathrm{rs}}=10z_{\mathrm{0}}$ ()). Here, we use $z_{*}=2.75h_{\mathrm{rs}}$.

Velocity scales (friction velocity, $u_{*}$, and convective velocity scale, $w_{*}$), Obukhov length ($L$), and planetary boundary layer height ($h$) characterising the stability regimes of LPDM-B simulations at measurement height $z_{\mathrm{m}}$ and with roughness length $z_{\mathrm{0}}$. Cases with measurement height within the roughness sublayer were disregarded (see text for details).

Scenario		$u_{*}$ [m s${}^{-\mathrm{1}}$]	$w_{*}$ [m s${}^{-\mathrm{1}}$]	$L$ [$\mathrm{m}$]	$h$ [$\mathrm{m}$]
1	convective	0.2	1.4	$-$15	2000
2	convective	0.2	1.0	$-$30	1500
3	convective	0.3	0.5	$-$650	1200
4	neutral	0.5	0.0	$\mathrm{\infty }$	1000
5	stable	0.4	–	1000	800
6	stable	0.4	–	560	500
7	stable	0.3	–	130	250
8	stable	0.3	–	84	200
Receptor heights at $z_{\mathrm{m}}/h=[0.005,0.01,0.075,0.25,0.50]$
Roughness lengths $z_{\mathrm{0}}=[0.01,0.1,0.3,1.0,3.0]$ m

As expected for such a broad range of scenarios, the resulting footprints of LPDM-B simulations show a vast range of extents and sizes. Figure  depicts this range by means of peak location of the footprints and their extent, when integrated from their peak to 80 % contribution of the total footprint (cf. Sect. ). For example, the 80 % footprint extent ranged from a few tens to a few hundreds of metres upwind of the tower location for the lowest measurement heights. For the highest measurements, the 80 % footprints ranged up to 270 km.

Range of peak locations ($x_{max}$) and extent of 80 % crosswind-integrated footprints ($x_{80}$) of LPDM-B simulations as in Table . The colour depicts the simulated measurement heights, symbols stand for modelled roughness lengths, $z_{\mathrm{0}}=0.01$ m ($\mathit{▵}$), 0.1 m ($\mathit{\square }$), 0.3 m ($\circ$), 1.0 m ($\mathrm{♢}$), 3.0 m ($\mathit{▿}$). Note the use of the log scale to accommodate the full range of the simulations.

[Figure omitted. See PDF]

An additional set of 27 LPDM-B simulations was run for independent evaluation of the footprint parameterisation. Measurement heights that are typical for flux tower sites were selected for this evaluation set, with boundary layer conditions again ranging from convective to stable. Table  lists the characteristics of these additional scenarios.

Scaling of footprints

The vast range in footprint sizes presented above clearly manifests that it is not practical to fit a single footprint parameterisation to all real-scale footprints. An additional step of footprint scaling is hence needed, with the goal of deriving a universal non-dimensional footprint. Ideally, such dimensionless footprints collapse to a single shape or narrow ensemble of curves. We follow a method that borrows from Buckingham $\mathrm{\Pi }$ dimensional analysis e.g., using dimensionless $\mathrm{\Pi }$ functions to scale the footprint estimates of LPDM-B, similar to , and scale the two components of the footprint, the crosswind-integrated footprint, and its crosswind dispersion (Eq. ), in two separate steps.

Scaled crosswind-integrated footprint

As in , we choose scaling parameters relevant for the crosswind-integrated footprint function, $\stackrel{\mathrm{‾}}{f^y}\left(x\right)$. The first choice is the receptor height, $z_{\mathrm{m}}$, as experience shows that the footprint (both its extent and footprint function value) is most strongly dependent on this height. Second, as indicated in Eq. (), the footprint is proportional to the flux at height $z_{\mathrm{m}}$. We hence formulate another scaling parameter based on the common finding that turbulent fluxes decline approximately linearly through the planetary boundary layer, from their surface value to the boundary layer height, $h$, where they disappear e.g.. Lastly, as a transfer function in turbulent boundary layer flow, the footprint is directly affected by the mean wind velocity at the measurement height, $\stackrel{\mathrm{‾}}{u}\left(z_{\mathrm{m}}\right)$, as well as by the surface shear stress, represented by the friction velocity, $u_{*}$. The well-known diabatic surface layer wind speed profile e.g. relates $\stackrel{\mathrm{‾}}{u}\left(z_{\mathrm{m}}\right)$ to the roughness length, $z_{\mathrm{0}}$, and the integrated form of the non-dimensional wind shear, ${\mathrm{\Psi }}_M$, that accounts for the effect of stability ($z_{\mathrm{m}}/L$) on the flow.

Velocity scales ($u_{*},w_{*}$), Obukhov length ($L$), planetary boundary layer height ($h$) describing the stability regimes of LPDM-B simulations for evaluation of the footprint parameterisation at measurement heights $z_{\mathrm{m}}$ and with roughness lengths $z_{\mathrm{0}}$.

Scenario		$u_{*}$ [m s${}^{-\mathrm{1}}$]	$w_{*}$ [m s${}^{-\mathrm{1}}$]	$L$ [$\mathrm{m}$]	$h$ [$\mathrm{m}$]
10	convective	0.20	1.00	$-$50	2500
11	convective	0.20	0.80	$-$100	2500
12	convective	0.25	0.75	$-$200	2160
15	neutral	0.40	0.00	$\mathrm{\infty }$	800
13	neutral	0.60	0.00	$\mathrm{\infty }$	1200
14	neutral	0.80	0.00	$\mathrm{\infty }$	1600
16	stable	0.30	–	200	310
17	stable	0.50	–	100	280
18	stable	0.35	–	50	170
$[z_{\mathrm{m}},z_{\mathrm{0}}]=$[20 m, 0.05 m], [30 m, 0.5 m], [50 m, 0.05 m]

With the above scaling parameters, we form four dimensionless $\mathrm{\Pi }$ groups as $$\begin{array}{cc}{\mathrm{\Pi }}_{\mathrm{1}}& =\stackrel{\mathrm{‾}}{f^y}{z}_{\mathrm{m}}\\ {\mathrm{\Pi }}_{\mathrm{2}}& =\frac{x}{z_{\mathrm{m}}}\\ {\mathrm{\Pi }}_{\mathrm{3}}& =\frac{h-z_{\mathrm{m}}}{h}=\mathrm{1}-\frac{z_{\mathrm{m}}}{h}\\ & {\mathrm{\Pi }}_{\mathrm{4}}& =\frac{\stackrel{\mathrm{‾}}{u}\left(z_{\mathrm{m}}\right)}{u_{*}}k=\ln \left(\frac{z_{\mathrm{m}}}{z_{\mathrm{0}}}\right)-{\mathrm{\Psi }}_M,\end{array}$$ where $k=0.4$ is the von Karman constant. We use ${\mathrm{\Psi }}_M$ as suggested by : $${\mathrm{\Psi }}_M=\left\{\begin{array}{ll}-5.3\frac{z_{\mathrm{m}}}{L}& \text{for}L>\mathrm{0},\\ \ln \left(\frac{\mathrm{1}+{\mathit{\chi }}^{\mathrm{2}}}{\mathrm{2}}\right)+\mathrm{2}\ln \left(\frac{\mathrm{1}+\mathit{\chi }}{\mathrm{2}}\right)\\ -\mathrm{2}{\tan }^{-\mathrm{1}}\left(\mathit{\chi }\right)+\frac{\mathit{\pi }}{\mathrm{2}}& \text{for}L<\mathrm{0}\end{array}\right.$$ with $\mathit{\chi }=(\mathrm{1}-19z_{\mathrm{m}}/L{)}^{\mathrm{1}/\mathrm{4}}$. In principle, ${\mathrm{\Psi }}_M$ is based on Monin–Obukhov similarity and valid within the surface layer. Hence, special care was taken in testing this scaling approach for measurements outside the surface layer (see below). In contrast to , the present study incorporates the roughness length directly in the scaling procedure: high surface roughness (i.e. large $z_{\mathrm{0}}$) enhances turbulence relative to the mean flow, and thus shortens the footprints. Here, $z_{\mathrm{0}}$ is either directly used as input parameter or is implicitly included through the fraction of $\stackrel{\mathrm{‾}}{u}\left(z_{\mathrm{m}}\right)/u_{*}$.

The non-dimensional form of the crosswind-integrated footprint, $F^{y*}$, can be written as a yet unknown function $\mathit{\phi }$ of the non-dimensional upwind distance, $X^{*}$. Thus $F^{y*}=\mathit{\phi }\left(X^{*}\right)$, with $X^{*}={\mathrm{\Pi }}_{\mathrm{2}}{\mathrm{\Pi }}_{\mathrm{3}}{\mathrm{\Pi }}_{\mathrm{4}}^{-\mathrm{1}}$ and $F^{y*}={\mathrm{\Pi }}_{\mathrm{1}}{\mathrm{\Pi }}_{\mathrm{3}}^{-\mathrm{1}}{\mathrm{\Pi }}_{\mathrm{4}}$ , such that $$\begin{array}{}& X^{*}& =\frac{x}{z_{\mathrm{m}}}\left(\mathrm{1},-,\frac{z_{\mathrm{m}}}{h}\right){\left(\frac{\stackrel{\mathrm{‾}}{u}\left(z_{\mathrm{m}}\right)}{u_{*}},k\right)}^{-\mathrm{1}},& & =\frac{x}{z_{\mathrm{m}}}\left(\mathrm{1},-,\frac{z_{\mathrm{m}}}{h}\right){\left(\ln ,\left(\frac{z_{\mathrm{m}}}{z_{\mathrm{0}}}\right),-,{\mathrm{\Psi }}_M\right)}^{-\mathrm{1}},& F^{y*}& =\stackrel{\mathrm{‾}}{f^y}{z}_{\mathrm{m}}{\left(\mathrm{1},-,\frac{z_{\mathrm{m}}}{h}\right)}^{-\mathrm{1}}\frac{\stackrel{\mathrm{‾}}{u}\left(z_{\mathrm{m}}\right)}{u_{*}}k,& & =\stackrel{\mathrm{‾}}{f^y}{z}_{\mathrm{m}}{\left(\mathrm{1},-,\frac{z_{\mathrm{m}}}{h}\right)}^{-\mathrm{1}}\left(\ln ,\left(\frac{z_{\mathrm{m}}}{z_{\mathrm{0}}}\right),-,{\mathrm{\Psi }}_M\right).\end{array}$$

As a next step, the above scaling procedure is applied to all footprints of Scenarios 1 to 8 (Table ) derived by LPDM-B. Despite the huge range of footprint extents (Fig. ), the resulting scaled footprints collapse into an ensemble of footprints of very similar shape, peak location, and extent (Fig. ). Hence, the new scaling procedure for crosswind-integrated flux footprints proves to be successful across the whole range of simulations, including the large range of surface roughness lengths and stability regimes.

Density plot of scaled crosswind-integrated footprints of LPDM-B simulations as in Table  (low density: light blue; high density, 100 times denser than low density: dark red). The footprint parameterisation ${\widehat{F}}^{y*}\left({\widehat{X}}^{*}\right)$ (cf. Eq. ) is plotted as black line.

[Figure omitted. See PDF]

Density plot of real-scale (left panel) and scaled (right panel) lateral dispersion of LPDM-B simulations as in Table ; low density: light blue, high density (100 times denser than low density): dark red. The parameterisation ${\widehat{\mathit{\sigma }}}_y^{*}\left({\widehat{X}}^{*}\right)$ (cf. Eq. ) is plotted as black line.

[Figure omitted. See PDF]

Scaled crosswind dispersion

Crosswind dispersion can be described by a Gaussian distribution function with ${\mathit{\sigma }}_y$ as the standard deviation of the crosswind distance e.g.. In contrast to vertical dispersion, Gaussian characteristics are valid for crosswind dispersion for the entire stability range and are even appropriate over complex surfaces e.g.. The three-dimensional particle dispersion of LPDM-B incorporates the Gaussian lateral dispersion cf. . Note that variations of the mean wind direction by the Ekman effect are neglected. Hence, assuming Gaussian characteristics for the crosswind dispersion function, $D_y$, in Eq. (), the flux footprint, $f(x,y)$, can be described as e.g. $$f(x,y)=\stackrel{\mathrm{‾}}{f^y}\left(x\right){\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi }}{\mathit{\sigma }}_y}}\exp \left(-,{\displaystyle \frac{y^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_y^{\mathrm{2}}}}\right).$$ Here, $y$ is the crosswind distance from the centreline (i.e. the $x$ axis) of the footprint. The standard deviation of the crosswind distance, ${\mathit{\sigma }}_y$, depends on boundary layer conditions and the upwind distance from the receptor.

Similar to the crosswind-integrated footprint, we aim to derive a scaling approach of the lateral footprint distribution. We choose $y$ and ${\mathit{\sigma }}_y$ as the relevant length scales, and combine them with two scaling velocities, the friction velocity, $u_{*}$, and the standard deviation of lateral velocity fluctuations, ${\mathit{\sigma }}_v$. In addition to the $\mathrm{\Pi }$ groups of Eq. (), we therefore set $$\begin{array}{cc}{\mathrm{\Pi }}_{\mathrm{5}}& =\frac{y}{z_{\mathrm{m}}}\\ {\mathrm{\Pi }}_{\mathrm{6}}& =\frac{{\mathit{\sigma }}_y}{z_{\mathrm{m}}}\\ & {\mathrm{\Pi }}_{\mathrm{7}}& =\frac{{\mathit{\sigma }}_v}{u_{*}}.\end{array}$$ In analogy to, for example, , we define a non-dimensional standard deviation of the crosswind distance, ${\mathit{\sigma }}_y^{*}$, proportional to ${\mathrm{\Pi }}_{\mathrm{6}}{\mathrm{\Pi }}_{\mathrm{7}}^{-\mathrm{1}}$. The non-dimensional crosswind distance from the receptor, $Y^{*}$, is linked to ${\mathit{\sigma }}_y^{*}$ through Eq. () and accordingly has to be proportional to ${\mathrm{\Pi }}_{\mathrm{5}}{\mathrm{\Pi }}_{\mathrm{7}}^{-\mathrm{1}}$. With that $$\begin{array}{}& Y^{*}& =p_{\mathrm{s}\mathrm{1}}\frac{y}{z_{\mathrm{m}}}\frac{u_{*}}{{\mathit{\sigma }}_v},& {\mathit{\sigma }}_y^{*}& =p_{\mathrm{s}\mathrm{1}}\frac{{\mathit{\sigma }}_y}{z_{\mathrm{m}}}\frac{u_{*}}{{\mathit{\sigma }}_v},\end{array}$$ where $p_{\mathrm{s}\mathrm{1}}$ is a proportionality factor depending on stability. Based on the LPDM-B results, we set $p_{\mathrm{s}\mathrm{1}}=min(\mathrm{1},{\left|z_{\mathrm{m}},/,L\right|}^{-\mathrm{1}}{10}^{-\mathrm{5}}+p)$, with $p=0.8$ for $L\le \mathrm{0}$ and $p=0.55$ for $L>\mathrm{0}$. In Fig.  unscaled distance from the receptor, $x$, and scaled, non-dimensional distance $X^{*}$ are plotted against the unscaled and scaled deviations of the crosswind distance, respectively (see Appendix  for information on the derivation of ${\mathit{\sigma }}_y$ from LPDM-B simulations). The scaling procedure is clearly successful, as the scaled deviation of the crosswind distance, ${\mathit{\sigma }}_y^{*}$, of all LPDM-B simulations collapse into a narrow ensemble when plotted against $X^{*}$ (Fig. , right-hand panel). For large $X^{*}$, the ensemble spread increases mainly due to increased scatter of LPDM-B simulations for distances far away from the receptor.

Flux footprint parameterisation

The successful scaling of both along-wind and crosswind shapes of the footprint into narrow ensembles within a non-dimensional framework provides the basis for fitting a parameterisation curve to the ensemble of scaled LPDM-B results. Like for the scaling approach, the footprint parameterisation is set up in two separate steps, the crosswind-integrated footprint, and its crosswind dispersion.

Crosswind-integrated footprint parameterisation

The ensemble of scaled crosswind-integrated footprints $F^{y*}\left(X^{*}\right)$ of LPDM-B is sufficiently coherent that it allows for fitting a single representative function to it. We choose the product of a power function and an exponential function as a fitting function for the parameterised crosswind-integrated footprint ${\widehat{F}}^{y*}\left({\widehat{X}}^{*}\right)$: $${\widehat{F}}^{y*}=a({\widehat{X}}^{*}-d{)}^b\exp \left({\displaystyle \frac{-c}{{\widehat{X}}^{*}-d}}\right).$$ Derivation of the fitting parameters, $a,b,c,d$, is dependent on the constraint that the integral of the footprint parameterisation (cf. Eq. ) must equal unity to satisfy the integral condition ${\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{F}^{y*}\left(X^{*}\right)\text{d}{X}^{*}=\mathrm{1}$ cf. . Hence, $$\begin{array}{cc}\underset{d}{\overset{\mathrm{\infty }}{\int }}{\widehat{F}}^{y*}\left({\widehat{X}}^{*}\right)\mathrm{d}{\widehat{X}}^{*}& =\mathrm{1}\\ & & =\underset{d}{\overset{\mathrm{\infty }}{\int }}a({\widehat{X}}^{*}-d{)}^b\exp \left(\frac{-c}{{\widehat{X}}^{*}-d}\right)\mathrm{d}{\widehat{X}}^{*}& & =a{c}^{b+\mathrm{1}}\mathrm{\Gamma }(-b-\mathrm{1}),\end{array}$$ where $\mathrm{\Gamma }\left(b\right)$ is the gamma function, and here $\mathrm{\Gamma }(-b-\mathrm{1})\equiv {\int }_{\mathrm{0}}^{\mathrm{\infty }}{t}^{-b-\mathrm{2}}\exp (-t)\mathrm{d}t$.

The footprint parameterisation (Eq. ) is fitted to the scaled footprint ensemble using an unconstrained nonlinear optimisation technique based on the Nelder–Mead simplex direct search algorithm . With that, we find $$\begin{array}{cc}a& =1.452\\ b& =-1.991\\ c& =1.462\\ & d& =\mathrm{0.136\hspace{0.17em}.}\end{array}$$

Figure  shows that the parameterisation of the crosswind-integrated footprint represents all scaled footprints very well. The goodness-of-fit of this single parameterisation to the ensemble of scaled footprints for all simulated measurement heights, stability conditions, and roughness lengths is evident from model performance metrics see, e.g., including the Pearson's correlation coefficient ($R$), the fractional bias (FB), the fraction of the parameterisation within a factor of 2 of the scaled footprints (FAC2), the geometric variance (VG), and the normalised mean square error (NMSE). Table  lists these performance metrics for the parameterisation of the full extent of the crosswind-integrated footprint curve, for the footprint peak location, and for the footprint peak value of the parameterisation against the corresponding scaled LPDM-B results. The fit can be improved even more, if the parameters are optimised to represent footprints of convective or neutral and stable conditions only (see Appendix ).

Performance of the footprint parameterisation evaluated against all scaled footprints of LPDM-B simulations as in Table . Performance of the crosswind-integrated footprint parameterisation was tested against the complete footprint curve ($F^{y*}$), the footprint peak location ($X_{max}^{*}$), and the footprint peak value ($F_{max}^{y*}$). The parameterisation for crosswind dispersion (${\widehat{\mathit{\sigma }}}_y^{*}$) was similarly tested against ${\mathit{\sigma }}_y^{*}\left(x\right)$ for all scaled footprints. In all cases the following performance measures were used: Pearson's correlation coefficient ($R$), fractional bias (FB), fraction of the parameterisation within a factor of 2 of the scaled footprints (FAC2), geometric variance (VG), and normalised mean square error (NMSE). Note that $R$ is not evaluated for $X_{max}^{*}$ and $F_{max}^{y*}$ as the footprint parameterisation provides a single value for each of these.

Performance metrics	$F^{y*}$	$X_{max}^{*}$	$F_{max}^{y*}$	${\mathit{\sigma }}_y^{*}$
$R$	0.96	–	–	0.90
FB	0.022	0.050	0.008	$-$0.120
FAC2	0.71	1.00	1.00	1.00
VG	2.86	1.04	1.05	1.02
NMSE	0.48	0.05	0.05	0.02

Parameterisation of the crosswind footprint extent

A single function can also be fitted to the scaled crosswind dispersion. In conformity with , the fitting function was chosen to be of the form $${\widehat{\mathit{\sigma }}}_y^{*}=a_{\mathrm{c}}{\left({\displaystyle \frac{b_{\mathrm{c}}({\widehat{X}}^{*}{)}^{\mathrm{2}}}{\mathrm{1}+c_{\mathrm{c}}{\widehat{X}}^{*}}}\right)}^{\mathrm{1}/\mathrm{2}}.$$ A fit to the data of scaled LPDM-B simulations results in $$\begin{array}{cc}{a}_{\mathrm{c}}& =2.17\\ b_{\mathrm{c}}& =1.66\\ & c_{\mathrm{c}}& =\mathrm{20.0\hspace{0.17em}.}\end{array}$$ The above parameterisation of the scaled deviation of the crosswind distance of the footprint is plotted in Fig.  (right panel). The performance metrics confirm that the ${\mathit{\sigma }}_y^{*}$ of the scaled LPDM-B simulations are very well reproduced by the parameterisation ${\widehat{\mathit{\sigma }}}_y^{*}$ (Table ).

Real-scale flux footprint

Typically, users of footprint models are interested in footprints given in a real-scale framework, such that distances (e.g. between the receptor and maximum contribution to the measured flux) are given in metres or kilometres. Depending on the availability of observed parameters, the conversion from the non-dimensional (parameterised) footprints to real-scale dimensions can be based on either Eqs. () and (), or on Eqs. () and (). For convenience, the necessary steps of the conversion are described in the following, by means of some examples.

Example footprint estimate for the convective Scenario 1 of Table , a measurement height of 20 m, and a roughness length of 0.01 m. The receptor is located at (0/0) m and the $x$ axis points towards the main wind direction. Footprint contour lines (left panel) are shown in steps of 10 % from 10 to 90 %.

[Figure omitted. See PDF]

Maximum footprint contribution

The distance between the receptor and the maximum contribution to the measured flux can be approximated by the peak location of the crosswind-integrated footprint. The maximum's position can be deduced from the derivative of Eq. () with respect to $X^{*}$: $${\widehat{X}}_{max}^{*}={\displaystyle \frac{-c}{b}}+d.$$ Using the fitting parameters as listed in Eq. () to evaluate ${\widehat{X}}_{max}^{*}$, the peak location is converted from the scaled to the real-scale framework applying Eq. () $$\begin{array}{cc}{x}_{max}& ={\widehat{X}}_{max}^{*}{z}_{\mathrm{m}}{\left(\mathrm{1},-,\frac{z_{\mathrm{m}}}{h}\right)}^{-\mathrm{1}}\frac{\stackrel{\mathrm{‾}}{u}\left(z_{\mathrm{m}}\right)}{u_{*}}k\\ & & =0.87z_{\mathrm{m}}{\left(\mathrm{1},-,\frac{z_{\mathrm{m}}}{h}\right)}^{-\mathrm{1}}\frac{\stackrel{\mathrm{‾}}{u}\left(z_{\mathrm{m}}\right)}{u_{*}}k,\end{array}$$ or alternatively applying Eq. () $$\begin{array}{cc}{x}_{max}& ={\widehat{X}}_{max}^{*}{z}_{\mathrm{m}}{\left(\mathrm{1},-,\frac{z_{\mathrm{m}}}{h}\right)}^{-\mathrm{1}}\left(\ln ,\left(\frac{z_{\mathrm{m}}}{z_{\mathrm{0}}}\right),-,{\mathrm{\Psi }}_M\right)\\ & & =0.87z_{\mathrm{m}}{\left(\mathrm{1},-,\frac{z_{\mathrm{m}}}{h}\right)}^{-\mathrm{1}}\left(\ln ,\left(\frac{z_{\mathrm{m}}}{z_{\mathrm{0}}}\right),-,{\mathrm{\Psi }}_M\right),\end{array}$$ with ${\mathrm{\Psi }}_M$ as given in Eq. (). Hence, $x_{max}$ can easily be derived from observations of $z_{\mathrm{m}}$, $h$, and $\stackrel{\mathrm{‾}}{u}\left(z_{\mathrm{m}}\right)$, $u_{*}$, or $z_{\mathrm{0}}$, $L$, and the constant value of ${\widehat{X}}_{max}^{*}$. For suggestions on how to estimate the planetary boundary layer height, $h$, if not measured, see Appendix .

Two-dimensional flux footprint

The two-dimensional footprint function can be calculated by applying the crosswind dispersion (Eq. ) to the crosswind-integrated footprint. With inputs of the scaling parameters $z_{\mathrm{m}}$, $h$, $u_{*}$, ${\mathit{\sigma }}_v$, and $\stackrel{\mathrm{‾}}{u}\left(z_{\mathrm{m}}\right)$ or $z_{\mathrm{0}}$, $L$, the two-dimensional footprint for any $(x,y)$ combination can be derived easily by the following steps:

• evaluate $X^{*}$ using Eqs. () or () for given $x$;

• derive ${\widehat{F}}^{y*}$ and ${\widehat{\mathit{\sigma }}}_y^{*}$ by inserting $X^{*}$ for ${\widehat{X}}^{*}$ in Eqs. () and ();

• invert Eqs. () or () and () to derive $\stackrel{\mathrm{‾}}{f^y}$ and ${\mathit{\sigma }}_y$, respectively;

• evaluate $f(x,y)$ for given $x$ and $y$ using Eq. ().

Figure  depicts an example footprint for convective conditions, computed by applying the described approach to arrays of $(x,y)$ combinations.

Relative contribution to the total footprint area

Often, the interest lies in the extent and location of the area contributing to, for example, 80 % of the measured flux. For such applications, there are two approaches: (i) the crosswind-integrated footprint function, $\stackrel{\mathrm{‾}}{f^y}\left(x\right)$, is integrated from the receptor location to the upwind distance where the contribution of interest is obtained; (ii) the two-dimensional footprint function $f(x,y)$ is integrated from the footprint peak location into all directions along constant levels of footprint values until the contribution of interest is obtained. The result is the source area: the smallest possible area containing a given relative flux contribution cf. . This approach can also be used as a one-dimensional equivalent to the source area, for the crosswind-integrated footprint.

For case (i) starting at the receptor location, we denote ${\widehat{X}}_R^{*}$ as the upper limit of the footprint parameterisation ${\widehat{F}}^{*}\left({\widehat{X}}^{*}\right)$ containing the area of interest, i.e. the fraction $R$ of the total footprint (that integrates to 1). The integral of Eq. () up to ${\widehat{X}}_R^{*}$ can be simplified (see Appendix  for details) as $$R=\exp \left({\displaystyle \frac{-c}{{\widehat{X}}_R^{*}-d}}\right).$$ With that, the distance between the receptor and ${\widehat{X}}_R^{*}$ can be determined very simply as $${\widehat{X}}_R^{*}={\displaystyle \frac{-c}{\ln \left(R\right)}}+d,$$ and in real scale $$\begin{array}{cc}& x_R=\left(\frac{-c}{\ln \left(R\right)},+,d\right)z_{\mathrm{m}}{\left(\mathrm{1},-,\frac{z_{\mathrm{m}}}{h}\right)}^{-\mathrm{1}}\frac{\stackrel{\mathrm{‾}}{u}\left(z_{\mathrm{m}}\right)}{u_{*}}k\\ & & =\left(\frac{-c}{\ln \left(R\right)},+,d\right)z_{\mathrm{m}}{\left(\mathrm{1},-,\frac{z_{\mathrm{m}}}{h}\right)}^{-\mathrm{1}}\left(\ln ,\left(\frac{z_{\mathrm{m}}}{z_{\mathrm{0}}}\right),-,{\mathrm{\Psi }}_M\right),\end{array}$$ where $R$ is a value between 0.1 and 0.9. As the above is based on the crosswind-integrated footprint, the derivation includes the full width of the footprint at any along-wind distance from the receptor.

There is no near-analytical solution for the description of the source area, the extent of the fraction $R$, when integrating from the peak location e.g.. Instead, the size of the source area has to be derived through iterative search. For crosswind-integrated footprints, the downwind (${\widehat{X}}_{R\mathrm{d}}^{*}<{\widehat{X}}_{max}^{*}$) and upwind (${\widehat{X}}_{max}^{*}<{\widehat{X}}_{R\mathrm{u}}^{*}$) distance from the receptor including the fraction $R$ can be approximated as a function of ${\widehat{X}}_R^{*}$, using LPDM-B results: $${\widehat{X}}_{R\mathrm{d},\mathrm{u}}^{*}=n_{\mathrm{1}}{\left({\widehat{X}}_R^{*}\right)}^{n_{\mathrm{2}}}+n_{\mathrm{3}}.$$ For the downwind limit ${\widehat{X}}_{R\mathrm{d}}^{*}$, $n_{\mathrm{1}}=0.44$, $n_{\mathrm{2}}=-0.77$, and $n_{\mathrm{3}}=0.24$. For the upwind limit, ${\widehat{X}}_{R\mathrm{u}}^{*}$, the approximation is split into two parts, $n_{\mathrm{1}}=0.60$, $n_{\mathrm{2}}=1.32$, and $n_{\mathrm{3}}=0.61$ for ${\widehat{X}}_{max}^{*}<{\widehat{X}}_R^{*}\le 1.5$, and $n_{\mathrm{1}}=0.96$, $n_{\mathrm{2}}=1.01$, and $n_{\mathrm{3}}=0.19$ for $1.5<{\widehat{X}}_R^{*}<\mathrm{\infty }$. The scaled distances ${\widehat{X}}_{*,R\mathrm{d}}$ and ${\widehat{X}}_{*,R\mathrm{u}}$ can again be transformed into real-scale values using Eqs. () or ().

If the size and position of the two-dimensional $R$-source area are of interest, but not the footprint function value (i.e. footprint weight) itself, the pairs of $x_R$ and $y_R$ describing its shape can be drawn from a lookup table of the scaled corresponding $X_R^{*}$ and $Y_R^{*}$ values. If the footprint function values are needed for weighting of source emissions or sinks, iterative search procedures have to be applied to each footprint. Figure  (left panel) illustrates examples of contour lines of $R$ fractions from 10 to 90 % of a footprint.

Footprint estimates for extended time series

The presented footprint model is computationally inexpensive and hence can be run easily for several years of data in, for example, half-hourly time steps. Each single data point can be associated with its source area by converting the footprint coordinate system to geographical coordinates, and positioning a discretised spatial array containing the footprint function onto a map or aerial image surrounding the receptor position. In many cases, an aggregated footprint, a so-called footprint climatology, is of more interest to the user than a series of footprint estimates. The aggregated footprint can be normalised and presented for several levels of relative contribution to the total aggregated footprint. Figure  shows an example of such a footprint climatology for 1 month of half-hourly input data for the ICOS flux tower site Norunda in Sweden cf. . A footprint climatology can be derived for selected hours of several days, for months, seasons, years, etc., depending on interest.

Example footprint climatology for the ICOS flux tower Norunda, Sweden, for 1–31 May 2011. The red dot depicts the tower location with a receptor mounted at $z_{\mathrm{m}}=z_{\mathrm{receptor}}-z_{\mathrm{d}}=12$ m. Footprint contour lines are shown in steps of 10 % from 10 to 90 %. The background map is tree height derived from an airborne lidar survey.

[Figure omitted. See PDF]

Combined with remotely sensed data, a footprint climatology provides spatially explicit information on vegetation structure, topography, and possible source/sink influences on the measured fluxes. This additional information has proven to be beneficial for analysis and interpretation of flux data e.g.. A combination of the footprint parameterisation presented here with high-resolution remote sensing data can be used not only to estimate the footprint area for measurements, but also to weigh or classify spatially continuous information on the surface and vegetation for its impact on measurements.

Certain remotely sensed data, for example airborne lidar data, allow for approximate derivation of the zero-plane displacement height and the surface roughness length . Alternatively, $z_{\mathrm{d}}$ and $z_{\mathrm{0}}$ may be estimated from flux tower measurements e.g.. If these measures vary substantially for different wind directions, we suggest running a spin-up of the footprint model, updating the measurement height ($z_{\mathrm{m}}=z_{\mathrm{receptor}}-z_{\mathrm{d}}$) and the footprint-weighted $z_{\mathrm{0}}$ input with each step, until a “steady state” of the footprints is reached. We recommend such a spin-up procedure despite the fact that footprint models are in principle not valid for non-scalars, such as momentum.

Performance of the footprint parameterisation evaluated against the second set of LPDM-B footprints of Table  (nine simulations for each stability regime), in real scale. Performance of the crosswind-integrated footprint parameterisation was tested against the crosswind-integrated footprint curve ($\stackrel{\mathrm{‾}}{f^y}$), the footprint peak location ($x_{max}$), the footprint peak value (${\stackrel{\mathrm{‾}}{f^y}}_{max}$), and the standard deviation of the crosswind distance at the peak location (${\mathit{\sigma }}_y\left(x_{max}\right)$). See Table  for abbreviations of performance measures.

Performance	$\stackrel{\mathrm{‾}}{f^y}$	$x_{max}$	${\stackrel{\mathrm{‾}}{f^y}}_{max}$	${\mathit{\sigma }}_y\left(x_{max}\right)$
metrics	[${\mathrm{m}}^{-\mathrm{1}}$]	[$\mathrm{m}$]	[${\mathrm{m}}^{-\mathrm{1}}$]	[$\mathrm{m}$]
Convective Scenarios 10, 11, 12
$R$	0.97	0.99	0.97	0.99
FB	0.014	0.109	$-$0.063	0.125
FAC2	0.99	1.00	1.00	1.00
VG	43.36	1.06	1.02	1.07
NMSE	0.33	0.09	0.03	0.07
Neutral Scenarios 13, 14, 15
$R$	0.97	0.99	0.96	0.47
FB	0.020	0.174	$-$0.067	0.05
FAC2	0.99	1.00	1.00	1.00
VG	1.62	1.11	1.04	1.12
NMSE	0.27	0.22	0.03	0.12
Stable Scenarios 16, 17, 18
$R$	0.99	0.99	0.99	0.88
FB	0.020	0.092	0.000	$-$0.013
FAC2	0.89	1.00	1.00	1.00
VG	1.16	1.03	1.03	1.02
NMSE	0.09	0.13	0.01	0.03

Discussion

Evaluation of FFP and sensitivity to input parameters

Exhaustive evaluation of footprint models is still a difficult task, and, clearly, tracer-flux field experiments would be very helpful. We are aware that in reality such experiments are both challenging and expensive to run. However, the aim of the present study is not to present a new footprint model, but to provide a simple and easily accessible parameterisation or “shortcut” for the much more sophisticated, but highly resource intensive, Lagrangian stochastic particle dispersion footprint model LPDM-B of . For the current study, we hence restrict the assessment of the presented footprint parameterisation to an evaluation against an additional set of LPDM-B simulations. A description of these additional scenarios can be found in Table .

The capability of the footprint parameterisation to reproduce the real-scale footprint of LPDM-B simulations is tested by means of the full extent of the footprint, its peak location, peak value, and its crosswind dispersion. Performance metrics show that for all stability classes (convective, neutral, and stable scenarios), the footprint parameterisation is able to predict the footprints simulated by the much more sophisticated Lagrangian stochastic particle dispersion model very accurately (Table ).

Results shown here clearly demonstrate that our objective of providing a shortcut to LPDM-B has been achieved. The full model was tested successfully against wind tunnel data . Further, the dispersion core of LPDM-B was evaluated successfully against wind tunnel and water tank data, LES results, and a full-scale tracer experiment . These considerations lend confidence to the validity of LPDM-B and thus FFP. They suggest that, despite its simplicity, FFP is suitable for a wide range of real-world applications, and is fraught with much less restrictive assumptions and turbulence regime limitations than what most other footprint models are faced with. We have applied the new scaling approach to LPDM-B, but it is likely similarly applicable to other complex footprint models.

Sensitivity of footprint peak location ($x_{max}$), peak value (${\stackrel{\mathrm{‾}}{f^y}}_{max}$), and the standard deviation of the crosswind distance at the peak location (${\mathit{\sigma }}_y\left(x_{max}\right)$) of the footprint parameterisation FFP to changes of the input parameters $h$ (boundary layer height) and $z_{\mathrm{0}}$ (roughness length) by $\pm$5, $\pm$10, $\pm$20 % for all scenarios of Table , in real scale. Changes are denoted in % deviation from the footprint parameterisation for the original input values of Table .

Change in input		$\mathrm{\Delta }{x}_{max}$	$\mathrm{\Delta }{\stackrel{\mathrm{‾}}{f^y}}_{max}$	$\mathrm{\Delta }{\mathit{\sigma }}_y\left(x_{max}\right)$
[%]		[%]	[%]	[%]
Convective Scenarios 10, 11, 12
$h$	$\pm$5	$\mp$0.1	$\pm$0.1	0.0
$h$	$\pm$10	$\mp$0.1	$\pm$0.1	0.0
$h$	$\pm$20	$\mp$0.3	$\pm$0.3	0.0
$z_{\mathrm{0}}$	$\pm$5	$\mp$1.1	$\pm$1.1	0.0
$z_{\mathrm{0}}$	$\pm$10	$\mp$2.2	$\pm$2.2	0.0
$z_{\mathrm{0}}$	$\pm$20	$\mp$4.4	$\pm$4.3	0.0
Neutral Scenarios 13, 14, 15
$h$	$\pm$5	$\mp$0.2	$\pm$0.2	0.0
$h$	$\pm$10	$\mp$0.3	$\pm$0.3	0.0
$h$	$\pm$20	$\mp$0.7	$\pm$0.7	0.0
$z_{\mathrm{0}}$	$\pm$5	$\mp$0.9	$\pm$0.9	0.0
$z_{\mathrm{0}}$	$\pm$10	$\mp$1.9	$\pm$1.9	0.0
$z_{\mathrm{0}}$	$\pm$20	$\mp$3.8	$\pm$3.7	0.0
Stable Scenarios 16, 17, 18
$h$	$\pm$5	$\mp$0.9	$\pm$0.9	$\pm$0.1
$h$	$\pm$10	$\mp$1.8	$\pm$1.7	$\pm$0.2
$h$	$\pm$20	$\mp$3.7	$\pm$3.6	$\pm$0.4
$z_{\mathrm{0}}$	$\pm$5	$\mp$0.7	$\pm$0.7	0.0
$z_{\mathrm{0}}$	$\pm$10	$\mp$1.4	$\pm$1.4	0.0
$z_{\mathrm{0}}$	$\pm$20	$\mp$2.8	$\pm$2.8	0.0

For the calculation of footprints with FFP, the values of the input parameters $z_{\mathrm{m}}$, $\stackrel{\mathrm{‾}}{u}\left(z_{\mathrm{m}}\right)$, $u_{*}$, $L$, and ${\mathit{\sigma }}_v$ can be derived from measurements typically available from flux towers. Input values for the roughness length, $z_{\mathrm{0}}$, may be derived from turbulence measurements or estimated using the mean height of the roughness elements e.g.. In the case of not perfectly homogeneous surfaces, these $z_{\mathrm{0}}$ values may vary depending on wind direction (see also Sect. ). Measurements of the boundary layer height, $h$, are available only rarely, and the accuracy of estimates of $h$ may vary substantially. In the following, we hence evaluate the sensitivity of the footprint parameterisation on the input parameters $z_{\mathrm{0}}$ and $h$.

The sensitivity of the FFP derived footprint estimate to changes in $h$ and $z_{\mathrm{0}}$ by $\pm$5, $\pm$10, and $\pm$20 % is tested for all scenarios of Table . For all scenarios, even changes of 20 % in $h$ and $z_{\mathrm{0}}$ result in only minor shifts or size alterations of the footprint (Table ). As to be expected, a small variation in $h$ does hardly alter footprint estimates for stability regimes with large $h$, namely, convective and neutral regimes. This finding is rather convenient, as reliable estimates of $h$ are difficult to derive for convective stabilities (see Appendix ). For stable scenarios, the footprint peak location is shifted closer to the receptor for overestimated $h$ and shifted further from the receptor for underestimated $h$. For these cases, overestimated $h$ will also very slightly increase the width of the footprint as described by ${\mathit{\sigma }}_y$ and vice versa. The impact of variations in the roughness length is quite similar for all atmospheric conditions, slightly decreasing the footprint extent for overestimated $z_{\mathrm{0}}$. Changes of the roughness length do not directly impact ${\mathit{\sigma }}_y$ but the absolute value of the footprint $f(x,y)$ can still vary, as a result of imposed changes in $\stackrel{\mathrm{‾}}{f^y}\left(x\right)$ (cf. Eq. ).

Comparison of FFP simulations for scenarios of Table  with corresponding simulations of the model (HKC00) of . Plotted are the extent of 80 % crosswind-integrated footprints, $x_{80}$ (left panel), and the crosswind dispersion at the footprint peak location, ${\mathit{\sigma }}_y\left(x_{max}\right)$ (right panel). The colour depicts the stability regime of the simulation: ML – mixed layer, FC – free convection layer, SL – surface layer (c for convective, n for neutral, and s for stable), NL – neutral layer, ZS – z-less scaling, and LS – local scaling seefor more details. Symbols denote modelled roughness lengths, $z_{\mathrm{0}}=0.01$ m ($\mathit{▵}$), 0.1 m ($\mathit{\square }$), 0.3 m ($\circ$), 1.0 m ($\mathrm{♢}$), 3.0 m ($\mathit{▿}$), and dashed lines denote 2 : 1 and 1 : 2, respectively. Note the use of the log–log scale to accommodate the full range of the simulations.

[Figure omitted. See PDF]

Limitations of FFP

Since FFP is based on LPDM-B simulations, LPDM-B's application limits are also applicable to FFP. As for most footprint models, these include the requirements of stationarity and horizontal homogeneity of the flow over time periods that are typical for flux calculations (e.g. 30–60 min). If applied outside these restrictions, FFP will still provide footprint estimates, but their interpretation becomes difficult and unreliable. Similarly, LPDM-B does not include roughness sublayer dispersion near the ground, nor dispersion within the entrainment layer at the top of the convective boundary layer. Hence, we suggest limiting FFP simulations to measurement heights above the roughness sublayer and below the entrainment layer (e.g. for airborne flux measurements). The $\mathrm{\Pi }$ functions of the scaling procedure also set some limitations to the presented footprint parameterisation (see below). Further, the presented footprint parameterisation has been evaluated for the range of parameters of Table  and application outside this range should be considered with care. For calculations of source areas of fractions $R$ of the footprint, we suggest $R\le 0.9$ (note that the source area for $R=\mathrm{1}$ is infinite). In most cases, $R=0.8$ is sufficient to estimate the area of the main impact to the measurement.

The requirements and limits of FFP for the measurement height and stability mentioned above can be summarised as follows: $$\begin{array}{cc}20z_{\mathrm{0}}& <z_{\mathrm{m}}<h_{\mathrm{e}}\\ & -15.5& \le \frac{z_{\mathrm{m}}}{L},\end{array}$$ where 20 $z_{\mathrm{0}}$ is of the same order as the roughness sublayer height, $z_{*}$ (see Sect. ), and $h_{\mathrm{e}}$ is the height of the entrainment layer typically, $h_{\mathrm{e}}\approx 0.8h$, e.g.. Equation () is required by ${\mathrm{\Pi }}_{\mathrm{4}}$ for a measurement height just above $z_{*}$ and may be adjusted for different values of $z_{*}$. At the same time, Eq. () also restricts application of the footprint parameterisation for very large measurement heights in strongly convective situations. For such cases, scaled footprints of LPDM-B simulations are of slightly shorter extent than those of the parameterisation and also include small contributions to the footprint from downwind of the receptor location (see Fig. ). Including the convective velocity scale as a scaling parameter did not improve the scaling. To account for such conditions, we hence suggest FFP parameters specific to the strongly convective stability regime (see Appendix ).

Same as Fig.  but for results of FFP compared with corresponding simulations of the model (KM01) of .

[Figure omitted. See PDF]

Comparison with other footprint models

In the following, we compare footprints of three of the a most commonly used models with results of FFP: the parameterisation of with crosswind extension of , the model of , and the footprint parameterisation of , hereinafter denoted HKC00, KM01, and KRC04, respectively.

For the comparison, the three above models and FFP were run for all scenarios listed in Table . As mentioned earlier, these scenarios span stability regimes ranging from the mixed layer (ML), free convection layer (FC), surface layer (SL; here further differentiated into convective, c, neutral, n, and stable, s), the neutral layer (NL), z-less scaling layer (ZS), and finally local scaling layer (LS). The wide range of stability regimes means that, unlike FFP, HKC00 and KM01 were in some cases run clearly outside their validity range. As, in practice, footprint models are run outside of their validity range quite frequently when they are applied to real environmental data, these simulations are included here.

Figure  shows the upwind extents of 80 % of the crosswind-integrated footprint ($x_{80}$) of these simulations of HKC00 against the corresponding results of FFP. Clearly, HKC00's footprints for neutral and stable scenarios extend further from the receptor than corresponding FFP's footprints by a factor of 1.5 to 2. This is the case for scenarios outside and within the surface layer. The results show most similar footprint extents for the convective part of the surface layer regime. In contrast, HKC00's footprints for elevated measurement heights within the free convection layer and footprints within the mixed layer are of shorter extent. The peak locations of the footprints (not shown) exhibit very similar behaviour to their 80 % extent. HKC00's crosswind dispersion is represented by ${\mathit{\sigma }}_y\left(x_{max}\right)$ at the footprint peak location (Fig. , right panel). HKC00's estimates are again larger than that of FFP for most scenarios except for mixed layer and free convection conditions; this is also found at half and at twice the peak location (not shown).

The along-wind extents of the footprint predictions of KM01 are very similar to HKC00's results, and hence the comparison of KM01 against FFP is similar as well: larger footprint extents resulting from KM01 than from FFP in most cases except for free convection and mixed layer scenarios, where FFP's footprints extend further (Fig. ). Again, the peak location of the footprints also follow this pattern (not shown). have discussed possible reasons for differences between KM01 and KRC04, relating these to LPDM-B capabilities of modelling along-wind dispersion that is also included in KRC04, but not in KM01. These reasons also apply to FFP. For crosswind dispersion, results of KM01 and FFP are relatively similar at the peak location of the footprint, $x_{max}$ (Fig. ). Differences are most evident for the neutral surface layer and for measurement heights above the surface layer. Nevertheless, the shape of the two-dimensional footprint is different between the two models. For most scenarios, the footprint is predicted to be wider by KM01 downwind of the footprint peak, and for scenarios within SL${}_{\text{c}}$ and FC it is predicted to be narrower upwind of the peak (not shown).

KRC04 and FFP were both developed on the basis of LPDM-B simulations. Hence, as expected, the results of these two footprint parameterisations agree quite well (Fig. ). FFP suggests that footprints extend slightly further from the receptor than KRC04 does, the difference is increasing with measurement height. FFP and KRC04 footprint predictions clearly differ for elevated measurement heights within the neutral layer, local scaling, and z-less scaling scenarios, which is due to the improved scaling approach of FFP. The footprint peak locations are predicted to be further away from the receptor by KRC04 than by FFP, with the difference decreasing for increasing measurement height (not shown).

To date, the availability of observational data suitable for direct evaluation of footprint models is very limited, and hence the performance of footprint models cannot be tested against “the truth”. Nevertheless, as stated in Sect. , LPDM-B and its dispersion core, the basis for FFP, have been evaluated successfully against experimental data, supporting the validity of FFP results.

Summary

Flux footprint models describe the area of influence of a turbulent flux measurement. They are typically used for the design of flux tower sites, and for the interpretation of flux measurements. Over the last decades, large monitoring networks of flux tower sites have been set up to study greenhouse gas exchanges between the vegetated surface and the lower atmosphere. These networks have created a great demand for footprint modelling of long-term data sets. However, to date available footprint models are either too slow to process such large data sets, or are based on too restrictive assumptions to be valid for many real-case conditions (e.g. large measurement heights or turbulence conditions outside Monin–Obukhov scaling).

In this study, we present a novel scaling approach for real-scale two-dimensional footprint data from complex models. The approach was applied to results of the backward Lagrangian stochastic particle dispersion model LPDM-B. This model is one of only few that have been tested against wind tunnel experimental data. LPDM-B's dispersion core was specifically designed to include the range from convective to stable conditions and was evaluated successfully using wind tunnel and water tank data, large-eddy simulation and a field tracer experiment.

The scaling approach forms the basis for the two-dimensional flux footprint parameterisation FFP, as a simple and accessible shortcut to the complex model. FFP can reproduce simulations of LPDM-B for a wide range of boundary layer conditions from convective to stable, for surfaces from very smooth to very rough, and for measurement heights from very close to the ground to high up in the boundary layer. Unlike any other current fast footprint model, FFP is hence applicable for daytime and night-time measurements, for measurements throughout the year, and for measurements from small towers over grassland to tall towers over mature forests, and even for airborne surveys.

Comparison of FFP simulations for scenarios of Table  with corresponding simulations of the footprint parameterisation (KRC04) of . Plotted is the extent of 80 % crosswind-integrated footprints ($x_{80}$). The colour depicts the stability regime of the simulation and symbols denote modelled roughness lengths (see Fig.  for details).

[Figure omitted. See PDF]

Footprint parameterisation optimised for specific stability conditions

There may be situations where footprint estimates are needed for only one specific stability regime, for example, when footprints are calculated for only a short period of time, or for a certain daytime over several days. For such cases, it may be beneficial to use footprint parameterisation settings optimised for this stability regime only. While scaled footprint estimates for neutral and stable conditions collapse to a very narrow ensemble of curves, footprints for strongly convective situations may also include contributions from downwind of the receptor location. For neutral and stable conditions, a specific set of fitting parameters for FFP has been derived using the LPDM-B simulations of Scenarios 4 to 6 (Table ). For convective conditions, additional LPDM-B simulations to Table  (Scenarios 1 to 3) have been included, to represent more strongly convective situations. These simulations (Scenario 1*) were run for the same set of receptor heights and surface roughness length as listed in Table , but with $u_{*}=0.2$ m s${}^{-\mathrm{1}}$, $w_{*}=2.0$ m s${}^{-\mathrm{1}}$, $L=-\mathrm{5}$ m, and $h=2000$ m. The resulting values for the parameters $a$, $b$, $c$, and $d$ for the crosswind-integrated footprint parameterisation, and $a_{\mathrm{c}}$, $b_{\mathrm{c}}$, and $c_{\mathrm{c}}$ for the crosswind dispersion parameterisation are listed in Table .

Please note that when applying these fitting parameters, the footprint functions for convective and neutral/stable conditions will not be continuous. We hence suggest to use the universal fitting parameters of Table A1 for cases where a transition between stability regimes may occur.

Fitting parameters of the crosswind-integrated footprint parameterisation and of the crosswind footprint extent for a “universal” regime (Scenarios 1 to 9), and for specifically convective (Scenarios 1* and 1 to 3) or neutral and stable regimes (Scenarios 4 to 9). For each scenario, all measurement heights and roughness lengths were included. See Table  and Appendix  for a description of the scenarios.

Stability	Universal	Convective	Neutral
regime			and stable
$a$	1.452	2.930	1.472
$b$	$-$1.991	$-$2.285	$-$1.996
$c$	1.462	2.127	1.480
$d$	0.136	$-$0.107	0.169
$a_{\mathrm{c}}$	2.17	2.11	2.22
$b_{\mathrm{c}}$	1.66	1.59	1.70
$c_{\mathrm{c}}$	20.0	20.0	20.0

Derivation of the boundary layer height

Determination of the boundary layer height, $h$, is a delicate matter, and no single “universal approach” can be proposed. Clearly, any available nearby observation (e.g. from lidar or radio sounding) should be used. $h$ may also be diagnosed from assimilation runs of high-resolution numerical weather predictions. For unstable (daytime) conditions, gave a comprehensive overview on different methods and discuss the associated caveats and uncertainties. For stable conditions, and provided a theoretical assessment of the boundary layer height under various limiting conditions. If none of the above measurements or approaches for $h$ are applicable, a so-called “meteorological pre-processor” may be used. A non-exhaustive suggestion for the latter is provided in the following.

For stable and neutral conditions there are simple diagnostic relations with which the boundary layer height can be estimated. proposed an interpolation formula for neutral to stable conditions: $$h={\displaystyle \frac{L}{3.8}}\left[-,\mathrm{1},+,{\left(\mathrm{1},+,2.28,{\displaystyle \frac{u_{*}}{fL}}\right)}^{\mathrm{1}/\mathrm{2}}\right],$$ where $L$ is the Obukhov length ($L=-u_{*}^{\mathrm{3}}\stackrel{\mathrm{‾}}{\mathrm{\Theta }}/\left(kg\right(\stackrel{\mathrm{‾}}{w^{\prime }{\mathrm{\Theta }}^{\prime }}{)}_{\mathrm{0}})$), $\stackrel{\mathrm{‾}}{\mathrm{\Theta }}$ is the mean potential temperature, $k=0.4$ the von Karman constant, $g$ the acceleration due to gravity, $(\stackrel{\mathrm{‾}}{w^{\prime }{\mathrm{\Theta }}^{\prime }}{)}_{\mathrm{0}}$ the surface (kinematic) turbulent flux of sensible heat, and $f=\mathrm{2}\mathrm{\Omega }\sin \mathit{\varphi }$ is the Coriolis parameter ($\mathit{\varphi }$ being latitude and $\mathrm{\Omega }$ the angular velocity of the Earth's rotation). Equation () is widely used in pollutant dispersion modelling e.g. and has the desired property to have limiting values corresponding to theoretical expressions. For large $L$ (i.e. approaching near-neutral conditions), Eq. () tends towards $$h_{\text{n}}=c_{\text{n}}{\displaystyle \frac{u_{*}}{\left|f\right|}}.$$ In their meteorological pre-processor, recommended $c_{\text{n}}=0.3$ corresponding to . Strictly speaking, Eq. () is valid only as long as the stability of the free atmosphere is also close to neutral, i.e. $10<N|f{|}^{-\mathrm{1}}<70$, where $N$ is the Brunt–Väisälä frequency above the boundary layer defined as $N=(g/\stackrel{\mathrm{‾}}{\mathrm{\Theta }}\mathrm{d}\mathrm{\Theta }/\mathrm{d}z{)}^{\mathrm{1}/\mathrm{2}}$ . If this is not the case, but $L$ is still very large, using $c_{\text{nn}}=1.36$ the asymptotic limit becomes $$h_{\text{nn}}=c_{\text{nn}}{\displaystyle \frac{u_{*}}{|fN{|}^{\mathrm{1}/\mathrm{2}}}}.$$ For strongly stable conditions, Eq. () approaches $$h_{\text{s}}=c_{\mathrm{s}}{\left({\displaystyle \frac{u_{*}L}{\left|f\right|}}\right)}^{\mathrm{1}/\mathrm{2}},$$ as already proposed by , who later termed this equation an “intermediate asymptote” for the more complete asymptotic expression derived for the stable boundary layer height by . They assign Eq. () to an applicability range of $\mathrm{4}\ll u_{*}L/\left|f\right|\ll 100$. The numerical parameters in Eq. () suggest $c_{\mathrm{s}}\approx$ 0.4 for the limit of very strong stability. suggested $c_{\text{s}}\approx 0.63$ and later proposed $c_{\text{s}}\approx 0.3$ based on LES results (here, both values are adapted for the above definition of $L$).

For convective conditions, the boundary layer height cannot be diagnosed due to the nearly symmetric diurnal cycle of the surface heat flux. It must therefore be integrated employing a prognostic expression, and starting at sunrise (before the surface heat flux first becomes positive) when the initial height is diagnosed using one of the above expressions. The slab model of is based on a simplified TKE (turbulence kinetic energy) equation including thermal and mechanical energy. The resulting rate of change for the boundary layer height is implicit in $h$ and may be solved iteratively or using $h\left(t_i\right)$ to determine the rate of change to yield $h\left(t_{i+\mathrm{1}}\right)$: $$\begin{array}{cc}\frac{\mathrm{d}h}{\mathrm{d}t}& =\frac{(\stackrel{\mathrm{‾}}{w^{\prime }{\mathrm{\Theta }}^{\prime }}{)}_{\mathrm{0}}}{\mathit{\gamma }}\left[\left(\frac{h^{\mathrm{2}}}{(\mathrm{1}+\mathrm{2}A)h-\mathrm{2}BkL}\right)\right.\\ & & {\left..+,\frac{C{u}_{*}^{\mathrm{2}}T}{\mathit{\gamma }g\left[\right(\mathrm{1}+A)h-BkL]}\right]}^{-\mathrm{1}}.\end{array}$$ Here, $\mathit{\gamma }$ is the gradient of potential temperature above the convective boundary layer. The latter is often not available for typical applications and can be approximated by a constant parameter (e.g. $\mathit{\gamma }\approx 0.01$ K m${}^{-\mathrm{1}}$, a typical mid-latitude value adopted from ). It must be noted, however, that $\mathrm{d}h/\mathrm{d}t$ in Eq. () is quite sensitive to this parameter. The model parameters, finally, $A=0.2,B=2.5$, and $C=\mathrm{8}$, are derived from similarity relations in the convective boundary layer that have been employed to find the growth rate of its height.

Addressing the finite nature of stochastic particle dispersion footprint models

For the parameterisation of the scaled footprints, using continuous functions, we need to address an issue that arises from the discrete nature of LPDM-B. Like in all stochastic particle dispersion models, the number of particles ($n$) released is necessarily finite. Despite a large $n$ (typically, $n\sim {10}^{\mathrm{6}}$ for simulations of this study), the distribution of particle locations is also finite, with a finite envelope. Hence, estimates of dispersion statistics become spurious near the particle envelope. In particular, the touchdown distribution statistics that form the basis of footprint calculations are truncated in close vicinity of the receptor, creating a “blind zone” of the footprint. The extent of this blind zone is related to the finite time, $T$, a particle takes to travel the vertical distance between source/sink and the receptor at $z_{\mathrm{m}}$. As the vertical dispersion scales with $u_{*}$, $T$ can be expressed as $T=p_{\mathrm{s}\mathrm{2}}{z}_{\mathrm{m}}/u_{*}$ where $p_{\mathrm{s}\mathrm{2}}$ is a proportionality factor depending on stability. In effect, no particle touchdown can be scored closer to the receptor than the horizontal travel distance for time $T$.

For the crosswind-integrated footprint, mean advection of the particle plume over time $T$ is the principal effect of the blind zone. This effect can be accounted for by a shift in $X^{*}$ by a constant distance, $d$, which is treated as a free parameter and determined by the fitting routine. For the crosswind dispersion of the footprint, the effect of the blind zone needs to be corrected for by the contribution to crosswind dispersion over time $T$, which is not accounted for in the source/sink particle touchdown distribution of LPDM-B. This contribution to the dispersion can be estimated as ${\mathit{\sigma }}_{y,\mathrm{0}}={\mathit{\sigma }}_v*T$, in accordance with Taylor's classical results for the near-source limit . Hence ${\mathit{\sigma }}_y$ of Eq. () becomes ${\mathit{\sigma }}_y={\mathit{\sigma }}_{y,\mathrm{0}}+{\mathit{\sigma }}_{y,\mathrm{res}}$ where the latter is the crosswind dispersion explicitly resolved by LPDM-B. Based on LPDM-B results, we set $p_{\mathrm{s}\mathrm{2}}=$ 0.35, 0.35, 0.5 for convective ($L<\mathrm{0}$), neutral ($L\to \mathrm{\infty }$), and stable conditions ($L>\mathrm{0}$), respectively.

Derivation of relative contribution to total footprint area

For the integration of the footprint parameterisation ${\widehat{F}}^{*}\left({\widehat{X}}^{*}\right)$ to an upper limit at ${\widehat{X}}_{*,R}$, we introduce the auxiliary variables $x_d={\widehat{X}}_{*}-d$ and $r={\widehat{X}}_{*,R}-d$. With that, and based on Eq. (), the integral of the footprint parameterisation up to $r$ can be expressed as $$R=\underset{d}{\overset{{\widehat{X}}_{*,R}}{\int }}{\widehat{F}}^{*}\left({\widehat{X}}^{*}\right)\mathrm{d}{\widehat{X}}^{*}=\underset{\mathrm{0}}{\overset{r}{\int }}a(x_d{)}^b\exp \left({\displaystyle \frac{-c}{x_d}}\right)\mathrm{d}{x}_d.$$ Substituting $c/x_d$ with $t$, the above integral can be solved as follows $$\begin{array}{cc}R& =\underset{\mathrm{\infty }}{\overset{c/r}{\int }}a{\left(\frac{c}{t}\right)}^b\exp (-t)\frac{-c}{t^{\mathrm{2}}}\mathrm{d}t\\ & & =a{c}^{b+\mathrm{1}}\underset{c/r}{\overset{\mathrm{\infty }}{\int }}{t}^{-b-\mathrm{2}}\exp (-t)\mathrm{d}t,\end{array}$$

$$\begin{array}{cc}R& =a{c}^{b+\mathrm{1}}\mathrm{\Gamma }\left(-,b,-,\mathrm{1},,,c,/,r\right)\\ & & =a{c}^{b+\mathrm{1}}\mathrm{\Gamma }\left(-,b,-,\mathrm{1},,,\frac{c}{{\widehat{X}}_{*,R}-d}\right).\end{array}$$ Here, $\mathrm{\Gamma }(-b-\mathrm{1},c/({\widehat{X}}_{*,R}-d\left)\right)$ is the upper incomplete gamma function defined as $\mathrm{\Gamma }(s,l)\equiv {\int }_l^{\mathrm{\infty }}{t}^{s-\mathrm{1}}\exp (-t)\mathrm{d}t$.As $b=-1.991\approx -\mathrm{2}$, $a\approx c$ (Eq. ), and $c^{b+\mathrm{1}}\approx c^{-\mathrm{1}}$, Eq. () can be further simplified to $$\begin{array}{cc}R& \approx \underset{c/r}{\overset{\mathrm{\infty }}{\int }}\exp (-t)\mathrm{d}t\\ & & =\exp \left(\frac{-c}{r}\right)=\exp \left(\frac{-c}{{\widehat{X}}_{*,R}-d}\right).\end{array}$$

Code availability

The code of the presented two-dimensional flux footprint parameterisation FFP and the crosswind-integrated footprint can be obtained from www.footprint.kljun.net. The code is available in several platform-independent programming languages. Check the same web page for an online version of the footprint parameterisation and for updates.

Acknowledgements

The authors would like to thank Anders Lindroth, Michal Heliasz, and Meelis Mölder for the Norunda flux tower data. We also thank three anonymous reviewers for their helpful suggestions. This research was supported by the UK Natural Environment Research Council (NERC, NE/G000360/1), the Commonwealth Scientific and Industrial Research Organisation Australia (CSIRO, OCE2012), by the LUCCI Linnaeus center of Lund University, funded by the Swedish Research Council, Vetenskapsrådet, the Austrian Research Community (OeFG), the Royal Society UK (IE110132), and by the German Helmholtz programme ATMO and the Helmholtz climate initiative for regional climate change research, REKLIM. Edited by: J. Kala