On Neglecting Free-Stream Turbulence in Numerical Simulation of the Wind-Induced Bias of Snow Gauges

Abstract

Numerical studies of the wind-induced bias of precipitation measurements assume that turbulence is generated by the interaction of the airflow with the gauge body, while steady and uniform free-stream conditions are imposed. However, wind is turbulent in nature due to the roughness of the site and the presence of obstacles, therefore precipitation gauges are immersed in a turbulent flow. Further to the turbulence generated by the flow-gauge interaction, we investigated the natural free-stream turbulence and its influence on precipitation measurement biases. Realistic turbulence intensity values at the gauge collector height were derived from 3D sonic anemometer measurements. Large Eddy Simulations of the turbulent flow around a chimney-shaped gauge were performed under uniform and turbulent free-stream conditions, using geometrical obstacles upstream of the gauge to provide the desired turbulence intensity. Catch ratios for dry snow particles were obtained using a Lagrangian particle tracking model, and the collection efficiency was calculated based on a suitable particle size distribution. The collection efficiency in turbulent conditions showed stronger undercatch at the investigated wind velocity and snowfall intensity below 10 mm h⁻¹, demonstrating that adjustment curves based on the simplifying assumption of uniform free-stream conditions do not accurately portray the wind-induced bias of snow measurements.

1. Introduction

Precipitation measurements obtained from catching type gauges are subject to systematic biases ascribable to both instrumental [1] and environmental sources [2]. The focus of this paper is on the catching bias induced by wind, the most important environmental error resulting in a lower amount of collected precipitation than what is expected in undisturbed conditions. Counting errors are instead related to the ability of the instrument to correctly measure the collected amount of precipitation [3] and are not addressed in this work.

Any object immersed in a wind field generates an aerodynamic response that depends on its geometry. Catching-type precipitation gauges themselves generate airflow deformations around their physical body, with significant acceleration and vertical velocity components arising especially above the gauge collector [4,5]. These airflow features have the potential to deviate the trajectories of the approaching hydrometeors, inducing a generally lower efficiency in the collection of precipitation than in the absence of wind. Snow measurements are particularly affected by the wind-induced bias because snowflakes are more sensitive than raindrops to the action of the airflow velocity components [6]. The wind effect on precipitation gauges is a well-known measurement issue, and is operationally mitigated either in the field, by employing wind shields [7], or in post-processing, by adjusting the measured data [8].

The traditional derivation of adjustment curves from field test experiments [9,10] has been accompanied by theoretical studies based on computational fluid dynamic (CFD) simulations [11,12]. Adjustment curves derived from field experiments have limitations due to their strict relationship with the local wind climatology, site morphology, type of precipitation (liquid, solid, and mixed), and accuracy of the reference configuration (the pit gauge for liquid and the double fence intercomparison reference for solid precipitation). On the other hand, numerical simulations are affected by the modelling assumptions adopted and by difficulties in reproducing the influencing variables and accounting for the natural features of both the precipitation process (particles microphysics) and the environmental conditions (wind variability and turbulence).

Among the relevant modelling assumptions, most literature studies based on numerical simulation [11,13,14] assume that turbulence is only generated by the interaction of the airflow with the gauge body, therefore imposing a steady and uniform incoming airflow, even though wind is known to be turbulent in nature. The roughness of the site and the presence of obstacles characterize the level of turbulence intensity at the gauge height, and, in operational conditions, precipitation gauges are immersed in a turbulent flow. Further to the local generation of turbulence due to the flow-gauge interaction, in this work we investigated the natural free-stream turbulence intensity and its influence on precipitation measurement biases. To the best of the author’s knowledge, only one attempt to assess the impact of free-stream turbulence on the wind-induced bias has been performed [15].

In that work, the role of free-stream turbulence in smoothing the airflow disturbance around a calyx-shaped precipitation gauge was investigated by means of CFD simulations and wind tunnel tests. Results revealed that both the normalized vertical and horizontal components of the flow velocity above the gauge collector were less accentuated in the turbulent free-stream configuration than in uniform free-stream conditions, due to the energy dissipation induced by turbulent fluctuations. However, in a CFD approach based on the unsteady Reynolds-averaged Navier–Stokes (URANS) equations, the analysis was limited to the effect of turbulence on the airflow deformation around and above the gauge, while the impact on the hydrometeors trajectories was not quantified [15].

In the present work, we used a time dependent CFD approach, based on the large eddy simulation (LES) modelling framework, and we focused on the impact of the uniform flow assumption on the assessment of the gauge collection efficiency (CE). The LES approach (at the cost of a significantly increased computational burden) is indeed more suitable to represent the inherent free-stream turbulence of the incoming wind, because the turbulent scales are numerically solved—rather than modelled—down to a threshold scale. Additionally, we implemented in the modelling scheme the turbulence intensity levels obtained from field measurements, using a 3D sonic anemometer. Finally, we used a Lagrangian Particle Tracking (LPT) model to quantify the impact of the free-stream turbulence on the CE of the gauge, to demonstrate that the role of free-stream turbulence is modulated by the particle size distribution (PSD) of the specific precipitation event, and thus depends on the precipitation intensity, consistently with the latest literature developments [16,17].

2. Materials and Methods

2.1. Calculation of a Realistic Wind Turbulence Intensity

Various studies, mainly focused on the evaluation of the efficiency of wind turbines, report estimates of the free-stream turbulence intensity of natural wind [18,19], generally measured at high elevation (larger than 10 m) above the ground. Precipitation gauges are installed with their collector positioned quite close to the ground surface (usually at about 2 m for gauges supported by poles, and about 0.5 m for gauges directly installed on the ground surface). With the aim of obtaining suitable turbulence intensity values to characterize the wind near the ground surface, 3D sonic anemometer measurements (kindly provided by Environmental Measurements Ltd.—EML^©, Newcastle, UK) from the Nafferton (UK) field test site (see Figure 1a) were analyzed. Data were composed of 38 minutes of high-frequency (20 Hz) wind measurements (Figure 1b). The moving average (Figure 1b) was calculated by testing different dimensions of the moving window, N, and checking that the Reynolds average approach was satisfied [20]. This requires that the mean values of the turbulent fluctuations $\overline{u^{\prime }\left(t\right)}$, $\overline{{\mathrm{v}}^{\prime }\left(t\right)}$ and $\overline{w^{\prime }\left(t\right)}$ in all directions (x, y, z) must be null. A trial-and-error procedure was adopted: after choosing a tentative value for N, the wind measurements were divided in wind classes (defined by their central value U_ref from 1 to 5.5 m s⁻¹ with bin size of 0.5 m s⁻¹), based on their average magnitude. For each wind class, the mean of the turbulent fluctuations was calculated, and the N value associated with the minimum of the average fluctuations was chosen (resulting in N = 125). Finally, the relative turbulence intensity values ($I_u$, $I_v$ and $I_w$) for each wind class were calculated as follows:

$$I_u=\frac{\sqrt{\overline{u^{\prime 2}}}}{U_{ref}}, I_v=\frac{\sqrt{\overline{{\mathrm{v}}^{\prime 2}}}}{U_{ref}} \mathrm{and} I_w=\frac{\sqrt{\overline{w^{\prime 2}}}}{U_{ref}}$$

(1)

2.2. CFD Simulations

The airflow pattern around a Geonor T200B (hereinafter, Geonor) precipitation gauge was obtained from CFD simulations, under both uniform and turbulent free-stream conditions. The Geonor gauge is a traditional weighing gauge characterized by a “chimney” shape, with an orifice diameter D = 0.160 m and height equal to 0.740 m. To reproduce realistic turbulence conditions near the gauge, three geometrical obstacles (square pillars, 0.3 × 0.3 × 7 m) were introduced in the computational domain upstream of the gauge (Figure 2a).

The model of the gauge geometry and the three pillars used to generate turbulence were prepared in the Standard Triangulation Language (STL) format. The computational domain of dimensions 40 × 20 × 20 m was realized with the x-axis orientated along the stream-wise direction, the y-axis along the crosswise direction, and the z-axis along the vertical direction. The center of the gauge collector was located at (0, 0, −2), and the square pillars were positioned at x = −12 m and spaced by 0.3 m. They extended four meters above and three meters below the elevation of the gauge collector.

The three-dimensional spatial domain was discretized using an unstructured hybrid hexahedral/prismatic finite volume mesh. The computational mesh for the simulation under uniform free-stream conditions was composed of twelve million cells, while by introducing the three pillars, the number of cells reached fourteen million (Figure 2b), with a significant increase in the computational burden. The quality of the mesh was checked by using the geometry parameters of orthogonality, skewness, and aspect ratio [21]. In both cases, the max skewness was lower than 3, the max non-orthogonality was lower than 55, and the maximum aspect ratio was lower than 4.1.

A refinement region (Figure 3a) was realized around the gauge and stretched downwind to improve the numerical solution of the turbulent wake generated by the gauge. Around the gauge body, the computational mesh was discretized with gradually increasing refinements while approaching the gauge surface (Figure 3b). Additional refinement layers were realized close to the surface, and the minimum dimension of the cells nearest to the geometry was about 0.35 mm, corresponding to 1/10 of the collector rim thickness.

To obtain high quality CFD results, we performed LESs. With this numerical approach, the turbulent scales of motion were solved down to the smaller size of the computational mesh. Moreover, under turbulent free-stream conditions, the LES was able to fully reproduce the development of the turbulence wake downstream of the three obstacles before impacting on the gauge body.

The fluid, air, was modelled as a Newtonian incompressible fluid with kinematic viscosity υ_a = 1.2 × 10⁻⁵ m² s⁻¹ and density ρ_a = 1.3 kg m⁻³ at a reference environmental temperature T_a = 0 °C. For each configuration, at the inlet of the computational domain (y–z plane), the undisturbed wind speed, U_ref, was imposed parallel to the x-axis and it was maintained uniform and constant in time, while a null gradient condition was set for pressure. At the outlet y–z plane, opposite to the inlet, a null gradient condition for the velocity and the atmospheric pressure was imposed, while the lateral surfaces of the domain were set as symmetry planes. The ground and the gauge body were assumed as impermeable surfaces with a no-slip condition. In all computational cells, initial conditions were imposed equal to U_ref for the velocity and equal to zero for the relative pressure.

The LES approach uses a filtering operation to define a threshold (Equation (2)) between the directly solved turbulent scales and those that are modelled. In Equation (3), $u_i$ is the ith component of velocity, $\overline{u_i}$ denotes the filtered-scale velocity that is directly solved by the model, and $u_i^{\prime }$ is the subgrid-scale velocity.

$$\Delta ={\left(\Delta x\Delta y\Delta z\right)}^{1/3}$$

(2)

$$u_i^{\prime }=u_i-\overline{u_i}$$

(3)

Using this notation, the filtering operation can be made general through the convolution integral between $u_i$ and a filter function G:

$$\overline{u_i}\left(x,t\right)={\displaystyle \int \int \int }G\left(x-\zeta ;\Delta \right)u_i\left(\zeta ,t\right)d^3\zeta$$

(4)

The Navier–Stokes continuity and momentum equations assume the following form:

$$\frac{\partial \overline{u_i}}{\partial x_i}=0$$

(5)

$$\frac{\partial \overline{u_i}}{\partial t}+\overline{u_j}\frac{\partial \overline{u_i}}{\partial x_j}=-\frac{1}{\rho }\frac{\partial \overline{p}}{\partial x_i}+\upsilon \frac{{\partial }^2\overline{u_i}}{\partial x_j{}^2}-\frac{\partial {\tau }_{ij}}{\partial x_j}$$

(6)

The Smagorinsky sub-grid scale approach was adopted to model the sub-grid stresses $({\tau }_{ij}$) with a gradient-diffusion process in a similar way to a molecular motion, in the form:

$${\tau }_{ij}=2{\upsilon }_T{S}_{ij}$$

(7)

where the strain rate tensor ($S_{ij}$) and kinematic eddy viscosity (${\upsilon }_T$) are:

$$S_{ij}=\frac{1}{2}\left(\frac{\partial \overline{u_i}}{\partial x_j}+\frac{\partial \overline{u_j}}{\partial x_i}\right)$$

(8)

$${\upsilon }_T={\left(C_S\Delta \right)}^2\sqrt{S_{ij}{S}_{ij}}$$

(9)

$C_S$ is the Smagorinsky dimensionless empirical parameter, set equal to 0.2.

2.3. Particle Tracking and Collection Efficiency

Based on the airflow simulation, the improved LPT model proposed in ref [22] was used to model the trajectories of snow particles when approaching the gauge collector. Dry snow particles were chosen because they are the most sensitive to the airflow characteristics, and especially to the acceleration and the vertical and transversal velocity components that are affected by turbulence. The cross-section area, density, and volume of the simulated particles were set in compliance with the power law formulation reported in the work of ref [23]. The modelled particle diameters were set equal to d = 0.25, 0.5, 0.75 mm and from 1 to 8 mm, with bin size of 1 mm. The starting positions of the particles were set on a regular grid located upstream and above the gauge collector. Their terminal velocity was obtained by imposing the equilibrium between the acting forces associated with the gravity acceleration and the aerodynamic drag coefficient. The drag coefficient formulation provided by ref [24] for crystals was employed. For each simulated particle, the initial velocity components were set equal to the undisturbed wind speed in the horizontal direction; the null value along the transversal direction and the terminal velocity, calculated as a function of the particle diameter, in the vertical direction.

For each particle diameter, the number of particles, $n\left(d, U_{ref}\right)$, which fall inside the gauge collector in disturbed conditions due to the aerodynamic response of the gauge body was calculated from the LPT model. For the two investigated free-stream configurations and for each particle size, we obtained the catch ratio r, defined as the ratio between the number of particles captured by the gauge collector in disturbed airflow conditions, $n\left(d,U_{ref}\right)$, and the maximum number of particles, captured in undisturbed conditions, $n_{max}\left(d,U_{ref}\right)$.

The integral of the catch ratio on the full range of diameters investigated provides the numerical CE, in the form:

$$\mathrm{CE}\left(U_{ref}\right)=\frac{{{\displaystyle \int }}_0^d{\rho }_p V_p n\left(d,U_{ref}\right) N\left(d\right) \mathrm{d}d}{{{\displaystyle \int }}_0^d{\rho }_p V_p n_{max}\left(d,U_{ref}\right) N\left(d\right) \mathrm{d}d}$$

(10)

where ${\rho }_p$ and $V_p$ are the density and volume of the particles, respectively, $n_{max}\left(d\right)$ is the number of particles which should fall inside the gauge collector if the airflow field was not affected by the aerodynamic response of the gauge body, and $N\left(d\right)$ is the total number of particles per unit volume of air and per unit size interval having volume equal to the sphere of diameter comprised between d and d + dd.

The method requires the assumption of a suitable PSD, i.e., the distribution of $N\left(d\right)$ (usually expressed in mm⁻¹ m⁻³) as a function of the particle diameter d (mm), over the specified range of diameters. In this work, the PSD formulation proposed in ref [25] for dry snow was used, having the characteristic exponential function proposed by ref [26], and obtained by fitting experimental observations:

$$N\left(d\right)=N_0{e}^{-\Lambda d}$$

(11)

where $N_0$ mm⁻¹ m⁻³ is the intercept and Λ mm⁻¹ is the slope of the linear form of the curve in a semi-log plot.

For a widespread mid-latitude rain, a constant value for N₀ and a relationship for Λ as a function of the rainfall intensity (RI mm h⁻¹) was found, as follows [26]:

$$N_0=8000; \Lambda =4.1 R{I}^{-0.21}$$

(12)

For solid precipitation, based on field experiments, the two parameters of the distribution as a function of the snowfall intensity (SI [mm h⁻¹]) are expressed as follows [25]:

$$N_0=3800 S{I}^{-0.87};\Lambda =2.55 S{I}^{-0.48}$$

(13)

3. Results

In Figure 4, the mean value of the turbulent fluctuations (circles) and the associated standard deviation (bars) are plotted for each wind class, as derived from the analysis of high-frequency anemometer measurements by setting N = 125. On the right-hand axis, the sample size (histogram) for each wind class is reported. This result reveals that turbulent fluctuations exhibited average values very close to zero, except for the wind classes characterized by a small sample size (wind speed > 5 m s⁻¹). As expected, along the two horizontal directions (x and y) the relative turbulence intensity values ($I_u$, $I_v$) were very similar, while for the vertical direction (z) the relative turbulent intensity values $(I_w$) were lower (Figure 5). In all cases, the relative turbulence intensity values decreased with increasing wind speed and seemed to have an asymptotical constant limit, although data for high wind speed classes were too few to be conclusive; this behavior is in line with literature observations, albeit referring to higher elevation measurements [19].

To limit the computational burden, for both the investigated free-stream conditions, simulations were performed by imposing a single wind speed, equal to U_ref = 6 m s⁻¹, at the inlet surface (upwind y, z plane of the computational domain) as a boundary condition. When the airflow overtook the three obstacles, its mean velocity magnitude decelerated and assumed the value of 2.5 m s⁻¹, which became the new reference wind speed (U_ref) for the turbulent free-stream configuration. Thanks to the scalability (low Reynolds dependence) of the mean flow fields under uniform free-stream conditions [11], the CFD results were rescaled by assuming U_ref = 2.5 m s⁻¹, to make the results comparable.

For the free-stream turbulent configuration, the decay along the longitudinal direction (x/D) of the three relative turbulence intensity profiles at the gauge collector elevation and the streamwise symmetry axis is shown in Figure 6. The longitudinal distance (x/D) between the three obstacles and the gauge was calibrated to obtain, at the reference wind speed U_ref = 2.5 m s⁻¹, the desired level of turbulence as measured at the Nafferton (U.K.) field test site as $I_u$ = 0.25.

When the free-stream turbulent flow impacts on a gauge body, an accelerated zone (Figure 7) and updraft and downdraft components arise above the collector of the gauge. In this configuration, the gauge collector was totally immersed in the turbulent flow generated by the three obstacles, as shown in Figure 8, where the normalized magnitude (U_mag/U_ref) of the instantaneous flow velocity in the horizontal plane (x,y) at the gauge collector elevation is reported. Along the transversal direction (y/D), the flow field started to become uniform many diameters away from the collector. As shown in Figure 9, the transversal component of the flow velocity (U_y/U_ref) revealed the typical characteristics of vortex shedding which occurred with a regular frequency and with alternate opposite velocity directions.

In addition to the work presented in ref [15], where the free-stream turbulence effect was investigated using an unsteady Reynolds-averaged Navier–Stokes approach (URANS), and results were provided in terms of attenuation of the airflow updraft and acceleration above the collector of the gauge, in the present work the particle–fluid interaction (using a one-way coupled model) was also addressed, and the two free-stream turbulence conditions were compared in terms of catch ratios and collection efficiency.

The LPT model was run upon the mean airflow fields obtained from the LES for each free-stream turbulence condition. The calculated catch ratios for each particle diameter are listed in

Table 1. Catch ratios obtained for each solid particle size d at U_ref = 2.5 m s⁻¹ using the LPT model and based on the large eddy simulation (LES) mean airflow fields for the Geonor precipitation gauge under uniform and turbulent free-stream conditions.

d (mm)	0.25	0.5	0.75	1	2	3	4	5	6	7	8
Uniform	0.190	0.393	0.449	0.492	0.590	0.659	0.708	0.734	0.757	0.787	0.787
Turbulent	0.164	0.164	0.216	0.289	0.590	0.718	0.777	0.820	0.849	0.869	0.892

. A more pronounced undercatch emerged for small particles (less than 2 mm) under turbulent free-stream conditions with respect to the uniform case, while the opposite occurred for larger particles (d > 2 mm). This was due to the greater aptitude of the small size particles to follow the turbulent velocity fluctuations, while larger particles are more inertial. In particular, the free-stream turbulence had two main effects: it reduced the aerodynamic effect of the wind–gauge interaction, with lower velocity components near the gauge body (as demonstrated in ref [15]) and introduced velocity fluctuations in all directions.

When the particle size is small, particles are more sensitive to the turbulent fluctuations, and therefore catch ratios in uniform free-stream conditions are larger than in turbulent conditions. With increasing the particle size, and therefore its terminal velocity, particle trajectories are less sensitive to the turbulent fluctuations. Moreover, in turbulent free-stream conditions, they cross a less disturbed airflow field so that catch ratios become larger than in uniform free-stream conditions.

To obtain an overall evaluation of the free-stream turbulence impact on the collection performance of the gauge, the overall CE was computed using Equation 10 with the assumed PSD formulation as proposed in ref [25].

The impact on the CE is indeed modulated by the shape of the PSD. Given the dependence of the PSD on the snowfall intensity, the result is highly dependent on the specific precipitation event. The strong relationship of the CE with the SI was demonstrated by refs 16 and 17 for solid and liquid precipitation, respectively. Following the same approach, the impact of the free-stream turbulence intensity on the CE was calculated here for a set of reference snowfall intensity values from 0.5 to 50 mm h⁻¹.

The CE curves obtained for the uniform and turbulent free-stream conditions at U_ref = 2.5 m s⁻¹ are plotted as a function of SI in Figure 10, where an inversion of the sign of the relative difference between the two curves is evident at about 10 mm h⁻¹. Beneath that threshold, neglecting the free-stream turbulence in the simulation would lead to an overestimation of the CE, with a resulting underestimation of the wind-induced error of precipitation measurements. The opposite occurs beyond 10 mm h⁻¹, with an underestimation of the CE that leads to an overestimation of the wind-induced error. Note that the typical snowfall rates are well below that threshold, and the impact of neglecting the free-stream turbulence intensity may lead to an overestimation of the collection efficiency of about 30% at SI = 1 mm h⁻¹ and U_ref = 2.5 m s⁻¹.

The inversion in the sign of the relative difference between the CE calculated under uniform and turbulent free-stream conditions can be explained by considering the relative weight of particles smaller than 2 mm (the inversion threshold in the catch ratios reported in Table 1) in the PSD (Figure 11). Indeed, at SI = 1 mm h⁻¹, particles smaller than 2 mm account for about 99.7% of the total number of particles in the unit volume of the atmosphere, which results in about 82.8% of the equivalent water volume. At 10 mm h⁻¹, however, these values become equal to 77% of the total number of particles and 7% of the equivalent water volume.

4. Discussion

The obtained results confirm that neglecting the free-stream turbulence intensity in the numerical simulation of the airflow field leads to a biased assessment of the wind-induced undercatch of catching type precipitation gauges. Further work is needed to quantify this bias for a wide range of wind speed and gauge shape configurations, which does require an extensive numerical effort due to the significant computational burden of the LES approach. This is beyond the scope of the present work, whose aim is to demonstrate the impact of neglecting the free-stream turbulence in the assessment of the collection efficiency, and to explain the underlying physical reasons based on a joint fluid-dynamics and particle–fluid interaction numerical approach. The specific setup chosen to demonstrate this phenomenon made no special assumptions that could limit the validity of the results and their extension to similar configurations of the wind speed, turbulence intensity, PSD, and gauge shape, except for the numerical value of the assumed parameters.

The derived dependence of the impact of turbulence on the snowfall intensity is in agreement with results in the literature and highlights the basic role of the PSD in modulating the collection efficiency, depending on the microphysical characteristics of the specific precipitation event. Additional studies of the link between the mathematical form and parameters of the PSD and precipitation intensity in various meteorological conditions would possibly enhance our understanding of such a modulation effect, leading to more detailed CFD simulations of the exposure problem and quantification of the associated measurement biases. The operational benefit would reside in the derivation of adjustment curves (also called transfer functions in the literature) as a function of both wind velocity and precipitation intensity, which would better reproduce the results of field experiments.

Further consequences are expected on the design and implementation of field experiments, where the traditional setup for the comparison of measurements obtained from test gauges and reference configurations (shielded gauges) should be complemented by including additional measurements obtained from high-resolution three-dimensional anemometers and suitable disdrometers. In this sense, the development of precipitation measurement instruments based on non-catching principles (typically optical) would benefit the contemporary measurement of the PSD and precipitation intensity.