Comment on: “Corrections to the Mathematical Formulation of a Backwards Lagrangian Particle Dispersion Model” by Gibson and Sailor (2012: Boundary-Layer Meteorology 145, 399–406)

We discuss the results of Gibson and Sailor (Boundary-Layer Meteorol 145:399–406, 2012) who suggest several corrections to the mathematical formulation of the Lagrangian particle dispersion model of Rotach et al. (Q J R Meteorol Soc 122:367–389, 1996). While most of the suggested corrections had already been implemented in the 1990s, one suggested correction raises a valid point, but results in a violation of the well-mixed criterion. Here we improve their idea and test the impact on model results using a well-mixed test and a comparison with wind-tunnel experimental data. The new approach results in similar dispersion patterns as the original approach, while the approach suggested by Gibson and Sailor leads to erroneously reduced concentrations near the ground in convective and especially forced convective conditions.


Introduction
Based on pioneering work of Thomson (1987) and Luhar and Britter (1989), Rotach et al. (1996) developed a novel Lagrangian particle dispersion model that simulates dispersion in unstable, stable and neutral atmospheric conditions, whereas others are only valid for a single condition. As with most Lagrangian models, the model of Rotach et al. (1996) also fulfills the well-mixed criterion (Thomson 1987).
Later, Kljun et al. (2002) used the model as a "dispersion module" of LPDM-b, a Lagrangian particle dispersion footprint model that itself later formed the basis of the flux footprint parametrization (FFP) in one and two dimensions (Kljun et al. 2004a(Kljun et al. , 2015. The dispersion model was also adapted and evaluated for use over urban areas (Rotach 2001;Rotach et al. 2004;Stöckl 2015). Gibson and Sailor (2012) suggested several corrections to the mathematical foundations in Rotach et al. (1996). Since many subsequent studies are based on this model, a critical examination of these corrections seems necessary and is undertaken in the following. To avoid repetition, the reader is directed to Rotach et al. (1996), Gibson and Sailor (2012), or Stöckl (2015) for the theoretical formulation of the model. Only the relevant parts are explained here. The following uses the nomenclature of Gibson and Sailor (2012) with standard notation for velocity fluctuation components (u, v, w) and (co-)variances (e.g., σ 2 u = uu). (2012) 2.1 Gaussian Streamwise Turbulence Gibson and Sailor (2012) note that in Rotach et al. (1996) the description of the Gaussian longitudinal velocity variance (Rotach et al. 1996's Eq. 30) was missing the power of two at uw G . They state that

Corrections Suggested by Gibson and Sailor
and while this is correct, it describes a simple typographical error in the article text. The model code has been correct since at least 1998 and hence this correction will not be discussed further (see also the Editor's footnote in Gibson and Sailor 2012).

Convective Streamwise Probability Current
Next, Gibson and Sailor (2012) point out an issue with a constant in the formulation of the convective streamwise probability current ϕ C u in Rotach et al. (1996). They correctly derive in their Eq. 23 and compare it to Eq. 21 of Rotach et al. (1996), where ρ = uw G σ u,G σ w,G is the correlation coefficient between streamwise and vertical velocity fluctuations. The argument of their exponential function has a numerator of 1, while Rotach et al. (1996) incorrectly list a numerator of 2. However, this error had also been corrected in the model code in the 1990s, hence it will not be discussed here either.
Furthermore, Gibson and Sailor (2012) state that the error function's argument (abbreviated byû in Eq. 22 of Rotach et al. 1996) should have an additional (1 − ρ 2 ) term in the denominator. However, this is not correct, since in Rotach et al. (1996) where V is the velocity covariance matrix. Given that v is independent of u and w, the (1, 1) matrix element of the inverse covariance matrix, V −1 11 , can be written as Substituting V −1 11 in Eq. 3 leads tô which is identical to the term inside the error function in Eq. 2. Hence, even though the expression for ϕ C u of Gibson and Sailor (2012) is correct, the expression stated in Rotach et al. (1996) is correct too.

Solenoidal Probability Current
The third correction suggested by Gibson and Sailor (2012) requires background information. All models based on the Langevin equation require what Gibson and Sailor (2012) denote the probability currents ϕ i , even if they are not always explicitly named so (Rodean 1996). In the following, i stands for the directional component index (1 and 3 in the two-dimensional version and 1, 2 and 3 in the three-dimensional version). In the model of Rotach et al. (1996), where ϕ C i denotes the convective term of the model and ϕ G i the neutral/stable term, linked by a transition function F. A third term ϕ * i is required to ensure that ϕ i → 0 for |u| → ∞ (Thomson 1987). This third term has to be solenoidal in velocity space, because ϕ i is derived from where P tot is the total (joint) probability density function (pdf) of the particles' velocity fluctuations, which is assumed to be equal to the pdf of the Eulerian velocity fluctuations (Thomson 1987). A solenoidal ϕ * i does not affect Eq. 7 (for details see Rotach et al. 1996), as by definition It is not possible to uniquely define ϕ * i in multi-dimensional models (Thomson 1987), where variables in one dimension depend on those in others, as is the case in the model of Rotach et al. (1996). Any function that fulfills the criteria above (solenoidal, lim |u|→∞ ϕ i = 0) can be used as ϕ * i . This non-uniqueness is a well-known, but so far, unsolved problem (Thomson and Wilson 2012). Note that the addition of the third, lateral dimension (i = 2) in subsequent studies (de Haan and Rotach 1998;Kljun et al. 2002;Rotach et al. 2004;Stöckl 2015) does not affect any of this, because v is independent of u and w, and ϕ * v = 0 (de Haan and Rotach 1998).
With the above, we move on to the third correction suggested by Gibson and Sailor (2012), who point out unit inconsistencies in the formulation of ϕ * i in Rotach et al. (1996), Namely, the arguments of neither the error function in Eq. 9a nor the second exponential function in Eq. 9b are dimensionless. Additionally, ϕ * w has units of s −1 instead of the required m −1 . Gibson and Sailor (2012) solve this by introducing factors to the arguments of the corresponding functions, thereby changing Eq. 9 to become with β 1 = 1/(2σ u,G ) and β 2 = 2/β 1 , henceforth called the Gibson-Sailor correction (GSC). Alternatively, using β 1 = 1 s m −1 and β 2 = 1 m s −1 formally also solves the unit inconsistencies and does not require changing the model code (suggestion by one of us, M. W. Rotach, in Gibson and Sailor 2012). The GSC does solve the unit inconsistency in the earlier version. However, it violates the requirement of a solenoidal ϕ * i (Eq. 8) and is therefore incorrect. This violation can be resolved by changing the dimensionless constant 2 in the numerator of β 2 to 1 instead, henceforth called the corrected GSC (cGSC); β 1 remains unchanged.
The accordingly modified Eq. 10 is solenoidal, because Since ∂ϕ * u ∂u is identical to ∂ϕ * w ∂w with opposite sign, the requirement of Eq. 8 is fulfilled. If the numerator of β 2 equals 2, as Gibson and Sailor (2012) suggest, the factor in the denominator of ∂ϕ * w ∂w in Eq. 11b is 4 instead of 2 and Eq. 8 is violated. The other requirement of lim |u|→∞ ϕ i = 0 is fulfilled by all three versions of ϕ * i (when substituted into ϕ i , Eq. 6). Substituting ϕ * i in Eq. 6, using either the original version Eq. 9, the GSC version Eq. 10, or the cGSC version, and then taking the limit of ϕ i , where each velocity fluctuation component u i approaches ±∞ separately (18 limits in total) shows this quite readily, using Aw A − Bw B = 0 in ϕ C i (Luhar and Britter 1989). Details on the derivation are omitted here for brevity and because the factors β 1 and β 2 do not influence the limits of ϕ i .
In summary, there is no unique solution for ϕ * i , both the original version Eq. 9 and the cGSC version herein can be used in the model, even though the former has unit inconsistencies. Those do not influence the well-mixed state of the model and can be formally fixed by adding two parameters of value 1 with correcting units (Gibson and Sailor 2012) while not changing the model code. The GSC version, however, should not be used for reasons described above.

Impact of the ϕ * i Modifications on Dispersion
As described in Sects. 2.1 and 2.2, the first two corrections suggested by Gibson and Sailor (2012) have no impact on model results. To describe the impact of the GSC (cf. Sect. 2.3) on the model results, a well-mixed test (as in Duman et al. 2014) was undertaken. A large number of particles (10 6 ) were initially uniformly distributed in height, and the dispersion simulation was run for 2 h (simulated time) with a timestep of 0.1 s. At the end of the simulation, the heights of the particles were binned into 100 equal height-ranges, and the number of particles in each bin was normalized by the expected number of particles per bin, given a uniform distribution. To fulfill the well-mixed criterion of Thomson (1987), the normalized concentration (i.e., particle density) has to be unity for all heights. Due to the stochastic nature of the model, exact unity could only be achieved in the limit of an infinite number of particles, hence a level of noise is expected. Different stability scenarios were run with the relevant scaling variables summarized in Table 1. An additional scenario with stable stratification was also investigated but yielded the same result as the neutral case, so that it is not explicitly discussed here. The result of this well-mixed test is shown in Fig. 1.
Forced convective 0.88 2.08 −133.3 0.2 700 The first two are taken from Kljun et al. (2015), while the third is from wind-tunnel experiments (Fedorovich et al. 1996), scaled to the atmosphere in Kljun et al. (2004b). u * is the friction velocity, w * is the convective velocity scale, L is the Obukhov length, z 0 is the roughness length, and z i refers to the planetary boundary-layer depth Fig. 1 Results of testing the well-mixed criterion for three different scenarios (a)-(c). For their description see Table 1. Shown is the normalized concentration as a function of the non-dimensional height z/z i , where z i is the planetary boundary-layer depth. The three lines represent different model runs where ϕ * i follows three different formulations: original (from Rotach et al. 1996), GSC (modified according to Gibson and Sailor 2012), and cGSC (changes as suggested in Sect. 2.3) For neutral conditions (Fig. 1a), and similarly for stable stratification (not shown), the exact formulation of ϕ * i described above does not influence the well-mixed test-or even the simulation outcome-at all, indicated by the three almost identical curves. This does not come as a surprise, since the transition function F and consequently ∂ F/∂z is zero at all heights for these conditions, reducing ϕ * i to zero as well, because ϕ * i depend linearly on ∂ F/∂z (Eq. 9 and Eq. 10). This behaviour is not general, and other formulations of ϕ * i (not considered here) may very well influence the model in stable and neutral conditions. Gibson and Sailor (2012) report that "a stable atmosphere (L = 100 m) showed less than 5% difference in peak magnitude of the crosswind integrated flux footprint" (comparing their formulation to the original in Rotach et al. 1996). However, the difference should be zero, and their result is most likely caused by an insufficient number (5 × 10 4 ) of particles, which lead to a too low signal-to-noise ratio.
In the convective case (Fig. 1b), the results of the run with the original ϕ * i and the cGSC version still coincide and are approximately unity at all heights, but the GSC version deviates from unity near the ground, indicating a violation of the well-mixed criterion. The effect is noticeable for small heights z/z i , because ∂ F/∂z (and consequently ϕ * i ) is largest near the ground. In convective conditions, F is unity everywhere except near the ground, where mechanically produced turbulence results in a velocity distribution that, with decreasing z/zi, progressively approaches a Gaussian distribution (F → 0) of the vertical velocity, hence producing a profile of its derivative that is highest for small z/z i and tends towards zero with increasing height (Rotach et al. 1996). This effect is very visible when comparing Fig. 1b to c, where, for forced convection, the mean wind speed is higher and hence ∂ F/∂z becomes zero for larger z/z i , resulting in a larger effect of the incorrect ϕ * i . To demonstrate the effect of the GSC in a practical example, a comparison with the forced convection wind-tunnel studies of Fedorovich et al. (1996) is show in Fig. 2, similar to Kljun et al. (2004b), who already compared the model of Rotach et al. (1996) with the same wind-tunnel data. Displayed are vertical profiles of a dimensionless concentration (see Kljun et al. 2004b). Each panel shows the model results for increasing distance from the source, all taken at the center of the plume. In each panel the three resulting profiles corresponding to the three versions of ϕ * i are plotted (original, GSC, and cGSC). When the model employs the original ϕ * i and the cGSC version, the concentration profiles appear  Kljun et al. (2004b), based on wind-tunnel data from Fedorovich et al. (1996). Shown are vertical profiles of dimensionless concentration at increasing dimensionless distance X * from the source. Note the varying scaling of the horizontal axes between panels. Measurements ( ) by Fedorovich and Thäter (2002). Explanation of the different lines in Fig. 1 similar, while the concentrations using the GSC version are markedly lower near the ground. These characteristics increase with distance from the source, and imply that the vertical dispersion with the GSC transports particles erroneously higher, which was already visible in Fig. 1c. It is noted that the GSC version can, depending on the distance from the source, reproduce the measurements better (Fig. 2, middle panels) or worse (Fig. 2, first and last panel). This indicates that, despite pronounced differences between GSC and the other two simulations, these are not the major reason (deficiency) in the model in accounting for an optimal reproduction of the measured concentrations.
In conclusion, the impact of an incorrect formulation of ϕ * i can be pronounced in convective conditions. For the version proposed by Gibson and Sailor (2012), the influence near the ground is especially large, which is unfortunate, considering that the concentration near or at the ground is probably of greatest interest in many studies.