In the atmospheric surface layer, the wavelength of the peak in the vertical velocity spectrum $\Lambda_w$ decreases with increasing stable stratification and proximity to the surface and this dependence constrains our ability to perform high-Reynolds-number large-eddy simulation (LES). Near the ground, the LES filter cutoff $\Delta_f$ is comparable to or larger than $\Lambda_w$ and as a result the subfilter-scale (SFS) fluxes in LES are always significant and their contribution to the total flux grows with increasing stability.

We use the three-dimensional turbulence data collected during the Horizontal Array Turbulence Study (HATS) field program to construct SFS fluxes and variances that are modelled in LES codes. Detailed analysis of the measured SFS motions shows that the ratio $\Lambda_w/\Delta_f$ contains the essential information about stratification, vertical distance above the surface, and filter size, and this ratio allows us to connect measurements of SFS variables with LES applications. We find that the SFS fluxes and variances collapse reasonably well for atmospheric conditions and filter widths in the range $\Lambda_w/\Delta_f = [0.2,15]$. The SFS variances are anisotropic and the SFS energy is non-inertial, exhibiting a strong dependence on the stratification, large-scale shear, and proximity to the surface. SFS flux decomposition into modified-Leonard, cross-, and Reynolds terms illustrates that these terms are of comparable magnitude and scale content at large $\Lambda_w/\Delta_f$. As $\Lambda_w/\Delta_f \rightarrow 0$, the SFS flux approaches the-ensemble-average flux and is dominated by the Reynolds term. Backscatter of energy from the SFS motions to the resolved fields is small in the bulk of the surface layer, less than 20% for $\Lambda_w/\Delta_{f} < 2$.

A priori testing of typical SFS models using the HATS dataset shows that the turbulent kinetic energy and Smagorinsky model coefficients $C_k$ and $C_s$ depend on $\Lambda_w/\Delta_f$ and are smaller than theoretical estimates based on the assumption of a sharp spectral cutoff filter in the inertial range. $C_k$ and $C_s$ approach zero for small $\Lambda_w/\Delta_f$. Much higher correlations between measured and modelled SFS fluxes are obtained with a mixed SFS model that explicitly includes the modified-Leonard term. The eddy-viscosity model coefficients still retain a significant dependence on $\Lambda_w/\Delta_f$ with the mixed model. A dissipation model of the form $\epsilon = C_\epsilon E_s^{3/2}/\Delta_f$ is not universal across the range of $\Lambda_w/\Delta_f$ typical of atmospheric LES applications. The inclusion of a shear-stability-dependent length scale (Canuto & Cheng 1997) captures a large fraction of the variation in the eddy-viscosity and dissipation model coefficients.