The spectral velocity tensor for homogeneous boundary-layer turbulence

We have postulated a simple model for the spectral tensor _ij (k) of an anisotropic, but homogeneous turbulent velocity field. It is a simple generalization of the spectral tensor _inf ^{ij piso}(k) for isotropic turbulence and we show how in the limit of isotropy, _ij (k) becomes equal to _inf ^{ij piso}(k). Whereas _inf ^{ij piso}(k) is determined entirely by one scalar function of k = ¦k¦, namely the energy spectrum, we need three independent scalar functions of k to specify _ij (k). We show how it is possible by means of the three stream-wise velocity component spectra to determine the three scalar functions in _ij (k) by solving two uncoupled, ordinary linear differential equations of first and second order. The analytic form of the component spectra each has a set of three parameters: the variance and the integral length scale of the velocity component and a dimensionless parameter, which governs the curvature of the spectrum in the transition domain from the inertial subrange towards lower wave numbers. When the three sets of parameters are the same, the three spectra correspond to isotropic turbulence and they are all interrelated and related to the energy spectrum. We show how it is possible to obtain these spectral forms in the neutral surface layer and in the convective boundary layer from data reported in the literature. The spectral tensor is used to predict the lateral coherences for all three velocity components and these predictions are compared with coherences obtained in two experiments, one using three masts at a horizontally homogeneous site in Denmark and one employing two aircraft flying in formation over eastern Colorado. Comparison shows reasonable agreement although with considerable experimental scatter.The National Center for Atmospheric Research is sponsored by the National Science Foundation.     

Turbulence decay in stratified and homogeneous marine layers

A spectral approach is applied to shear-induced turbulence in stratified layers. A system of spectral equations for stationary balance of turbulent energy and temperature variances was deduced in the vicinity of the local shear scale L_U = (ε/U_Z³)^1/2. At wavenumbers between the inertial-convective (k^−5/3) and wak turbulence (k⁻³) subranges, additional narrow spectral intervals—‘production’ subranges—may appear (E k⁻¹, E_T k⁻²). The upper boundary of these subranges is determined as L_U, and the lower boundaries as L_R (ε/U_ZN²)^1/2(χ/T_Z²). It is shown that the scale L_U is a unique spectral scale that is uniform up to a constant value for every hydrophysical field. It appears that the spectral scale L_U is equivalent to the Thorpe scale L_Th for the active turbulence model. Therefore, if turbulent patches are generated in a background of permanent mean shear, a linear relation between temperature and mass diffusivities exists. In spectral terms, the fossil turbulence model corresponds to the regime of the Boldgiano-Obukhov buoyancy subrange (E k^−11/5, E_T k^−7/5). During decay the buoyancy subrange is expanded to lower and higher wavenumbers. At lower wavenumbers the buoyancy subrange is bounded by L_** = 3(χ^1/2/N^1/2T_Z), which is equivalent to the Thorpe scale L_Th. In such a transition regime only, when the viscous dissipation rate is removed from the set of main turbulence parameters, the Thorpe scale does not correlate with the buoyancy scale L_N ε^1/2/N^3/2 and fossil turbulence is realized. Oceanic turbulence measurements in the equatorial Pacific near Baker Island confirm the main ideas of the active and fossil turbulence models.     

On the similarity of the characteristics of turbulence in an unstable boundary layer

The statistics of turbulence, such as the standard deviation of fluctuating velocities, in an unstable atmospheric boundary layer are assumed to be characterized by the combination of three specific lengths, Monin-Obukhov length L, observation height z and the height of mixing layer h. Unlike Monin-Obukhov similarity, even near the ground the effect of h is taken into account. According to observation, the length scale of the vertical velocity is proportional to z at least near the ground, but the lateral component depends mostly on h alone. The length scale of the longitudinal component depends on z and h.     

Effects of Topographical Slope Angle and Atmospheric Stratification on Surface-Layer Turbulence

The effect of topographical slope angle and atmospheric stratification on turbulence intensities in the unstably stratified surface layer have been parameterized using observations obtained from a three-dimensional sonic anemometer installed at 8 m height above the ground at the Seoul National University (SNU) campus site in Korea for the years 1999–2001. Winds obtained from the sonic anemometer are analyzed according to the mean wind direction, since the topographical slope angle changes significantly along the azimuthal direction. The effects of the topographical slope angle and atmospheric stratification on surface-layer turbulence intensity are examined with these data. It is found that both the friction velocity and the variance for each component of wind normalized by the mean wind speed decrease with increase of the topographical slope angle, having a maximum decreasing rate at very unstable stratification. The decreasing rate of the normalized friction velocity (u _* /U) is found to be much larger than that of the turbulence intensity of each wind component due to the reduction of wind shear with increase in slope angle under unstable stratification. The decreasing rate of the w component of turbulence intensity (σ _w /U) is the smallest over the downslope surface whereas that of the u component (σ _u /U) has a minimum over the upslope surface. Consequently, σ _w /u _* has a maximum increasing rate with increase in slope angle for the downslope wind, whereas σ _u /u _* has its maximum for the upslope wind. The sloping terrain is found to reduce both the friction velocity and turbulence intensity compared with those on a flat surface. However, the reduction of the friction velocity over the sloping terrain is larger than that of the turbulence intensity, thereby enhancing the turbulence intensity normalized by the friction velocity over sloping terrain compared with that over a flat surface.     

The effect of idealized water waves on the turbulence structure and kinetic energy budgets in the overlying airflow

The influence of an idealized moving wavy surface on the overlying airflow is investigated using direct numerical simulations (DNS). In the present simulations, the bulk Reynolds number is Re = 8000 (; where U₀ is the forcing velocity of the flow, h the height of the domain and v the kinematic viscosity) and the phase speed of the imposed waves relative to the friction velocity, i.e., the wave age varies from very slow to fast waves. The wave signal is clearly present in the airflow up to at least 0.15λ (where λ is the wave length) and is present up to higher levels for faster waves. In the kinetic energy budgets, pressure transport is mainly of importance for slow waves. For fast waves, viscous transport and turbulent transport dominate near the surface. Kinetic energy budgets for the wave and turbulent perturbations show a non-negligible transport of turbulent kinetic energy directed from turbulence to the wave perturbation in the airflow. The wave-turbulent energy transport depends on the size, tilt, and phase of the wave-induced part of the turbulent Reynolds stresses.According to the DNS data, slow waves are more efficient in generating isotropic turbulence than fast waves.Despite the differences in wave-shape as well as in Reynolds number between the idealized direct numerical simulations and the atmosphere, there are intriguing similarities in the turbulence structure. Important information about the turbulence above waves in the atmosphere can be obtained from DNS—the data must, however, be interpreted with care.     

黄土高原复杂地形受中尺度运动影响的稳定边界层湍流特征

利用兰州大学半干旱气候与环境观测站(Semi-Arid Climate and Environment Observatory of Lanzhou University,简称SACOL)2008年12月观测资料,研究了稳定边界层湍流特征.使用涡动相关资料研究湍流通量时,定义湍流的平均时间τ内的中尺度运动是造成湍流统计量变化范围大的主要原因,稳定情形? τ取几十秒至几分钟.对梯度理查森数大于0.3的强稳定情形的湍流尺度分解(MRD)谱分析表明,感热通量在112.4～449.9 s存在谱隙,尺度大于谱隙的中尺度运动造成了通量观测资料离散性大,甚至有支配性影响.动量通量的谱隙在112.4～224.9 s之间.弱风时,中尺度运动的影响更大,垂直风速标准差以0.1的比率随中尺度风速变化;垂直风速标准差同广义风速表现出很好的相关性,并随着广义风速消失而消失.三维风速标准差与摩擦速度呈很好的线性关系,垂直、水平、横风风速的无量纲标准差分别为1.35、2.54、2.21.对湍流动能的研究发现,在梯度理查森数大于0.3的条件下,仍然存在连续的湍流.以湍动能为依据,分析了湍流的平稳时间长度,其长度随稳定度变化而变化,2008年12月7～11日从133.5 s变化到856.2 s,湍流平稳时间长度反映了中尺度运动的发生频率.     

Intermittency of Turbulence in the Atmospheric Boundary Layer: Scaling Exponents and Stratification Influence

We show the relationship between the intermittency of turbulence and the type of stratification for different atmospheric situations during the SABLES98 field campaign. With this objective, we first demonstrate the scaling behaviour of the velocity structure functions corresponding to these situations; next, we analyze the curvature of the scaling exponents of the velocity structure functions versus the order of these functions (ζ _p vs. p), where ζ _p are the exponents of the power relation for the velocity structure function with respect to the scale. It can be proved that this curve must be concave, under the assumption that the incompressible approximation does not break down at high Reynolds numbers. The physical significance of this kind of curvature is that the energy dissipation rate increases as the scale of the turbulent eddies diminishes (intermittency in the usual sense). However, the constraints imposed by stability, preventing full development of the turbulence, allow the function ζ _p versus p to show any type of curvature. In this case, waves of high frequency trapped by the stability, or bursts of turbulence caused by the breaking up of internal waves, may produce a redistribution of energy throughout the scaling range. Due to this redistribution, the variation with the scale of the energy dissipation rate may be smaller (decreasing the intermittency) and, even in more stable situations, this rate may diminish (instead of increasing) as the scale diminishes (convex form of the curve ζ _p vs. p).     

Experimental studies of the development of quasi-two-dimensional turbulence in stably stratified fluid

Quasi-two-dimensional turbulence was generated by towing an array of vertical cylinders through a tank which was filled with a two-layer stratified fluid. Sugar and Epsom salts were used, to give matching refractive indices for the two layers. The interface between the two layers was seeded with approximately 1000 neutrally buoyant particles. The evolution of this quasi-two-dimensional turbulence was visualized by photographing the fluorescent particles illuminated by a horizontal laser sheet traversing in the vertical direction. The three-dimensional particle velocity was obtained by digitizing the streaks. The evolution of the velocity correlations, length scales, one-dimensional and two-dimensional velocity and vorticity spectra were obtained for N = 5.72 s⁻¹, N = 4.43 s⁻¹, and N = 2.55 s⁻¹ (where N is the Brunt-Väisälä frequency). The results showed the physical process of inverse energy cascading and the formation of dominant vortical structures under the influence of density stratification. Compared with idealized two-dimensional turbulence, the flow is highly dissipative at high N, as a result of the frictional dissipation between the interface and the unstratified layers.     

Simplification and Validation of a Spectral-Tensor Model for Turbulence Including Atmospheric Stability

A spectral-tensor model of non-neutral, atmospheric-boundary-layer turbulence is evaluated using Eulerian statistics from single-point measurements of the wind speed and temperature at heights up to 100 m, assuming constant vertical gradients of mean wind speed and temperature. The model has been previously described in terms of the dissipation rate \(\epsilon \), the length scale of energy-containing eddies \(\mathcal {L}\), a turbulence anisotropy parameter \(\varGamma \), the Richardson number Ri, and the normalized rate of destruction of temperature variance \(\eta _\theta \equiv \epsilon _\theta /\epsilon \). Here, the latter two parameters are collapsed into a single atmospheric stability parameter z / L using Monin–Obukhov similarity theory, where z is the height above the Earth’s surface, and L is the Obukhov length corresponding to \(\{Ri,\eta _\theta \}\). Model outputs of the one-dimensional velocity spectra, as well as cospectra of the streamwise and/or vertical velocity components, and/or temperature, and cross-spectra for the spatial separation of all three velocity components and temperature, are compared with measurements. As a function of the four model parameters, spectra and cospectra are reproduced quite well, but horizontal temperature fluxes are slightly underestimated in stable conditions. In moderately unstable stratification, our model reproduces spectra only up to a scale \(\sim \) 1 km. The model also overestimates coherences for vertical separations, but is less severe in unstable than in stable cases.     

Direct Numerical Simulation of Stable Channel Flow at Large Stability

We consider a model for the stable atmospheric boundary at large stability, i.e. near the limit where turbulence is no longer able to survive. The model is a plane horizontally homogeneous channel flow, which is driven by a constant pressure gradient and which has a no-slip wall at the bottom and a free-slip wall at the top. At the lower wall a constant negative temperature flux is imposed. First, we consider a direct numerical simulation of the same channel flow. The simulation is computed with the neutral channel flow as initial condition and computed as a function of time for various values of the stability parameter h/L, where h is the channel height and L is related to the Obukhov length. We find that a turbulent solution is only possible for h/L < 1.25 and for larger values turbulence decays. Next, we consider a theoretical model for this channel flow based on a simple gradient transfer closure. The resulting equations allow an exact solution for the case of a stationary flow. The velocity profile for this solution is almost linear as a function of height in most of the channel. In the limit of infinite Reynolds number, the temperature profile has a logarithmic singularity at the upper wall of the channel. For the cases where a turbulent flow is maintained in the numerical simulation, we find that the velocity and temperature profiles are in good agreement with the results of the theoretical model when the effects of the surface layer on the exchange coefficients are taken into account. Frans Nieuwstadt, a recently retired member of the BLM Editorial Board and a well-known member of the boundary-layer/turbulence community, died unexpectedly on 18 May 2005. An obituary will appear in a later issue of BLM.     

Analysis of Turbulence Collapse in the Stably Stratified Surface Layer Using Direct Numerical Simulation

The nocturnal atmospheric boundary layer (ABL) poses several challenges to standard turbulence and dispersion models, since the stable stratification imposed by the radiative cooling of the ground modifies the flow turbulence in ways that are not yet completely understood. In the present work we perform direct numerical simulation of a turbulent open channel flow with a constant (cooling) heat flux imposed at the ground. This configuration provides a very simplified model for the surface layer at night. As a result of the ground cooling, the Reynolds stresses and the turbulent fluctuations near the ground re-adjust on times of the order of L/u _τ , where L is the Obukhov length scale and u _τ is the friction velocity. For relatively weak cooling turbulence survives, but when Re_L=Lu_t/n <~100{Re_L=Lu_\tau/\nu \lesssim 100} turbulence collapses, a situation that is also observed in the ABL. This criterion, which can be locally measured in the field, is justified in terms of the scale separation between the largest and smallest structures of the dynamic sublayer.     

Different Turbulent Regimes and Vertical Turbulence Structures of the Urban Nocturnal Stable Boundary Layer

Turbulence in the nocturnal boundary layer(NBL) is still not well characterized, especially over complex underlying surfaces. Herein, gradient tower data and eddy covariance data collected by the Beijing 325-m tower were used to better understand the differentiating characteristics of turbulence regimes and vertical turbulence structure of urban the NBL. As for heights above the urban canopy layer(UCL), the relationship between turbulence velocity scale(V_TKE) and wind speed(V) was con...     

Turbulence in an almond orchard: Spatial variations in spectra and coherence

The spatial variability of turbulence in a fully-leafed almond orchard was studied. Two three-dimensional sonic anemometers were used to measure turbulence spectra and coherence at different vertical and lateral separations inside the canopy. Peak frequencies of the horizontal velocity components, normalized by local horizontal wind speed, are greater in the canopy crown than in the trunkspace. Peak-normalized frequencies for the vertical velocity power spectra are similar in the canopy crown and in the subcanopy trunkspace. Spectral slopes in the inertial subrange are more negative than those predicted with Kolmogorov's -2/3 theory. It is thought that the foliage elements act to short-circuit the eddy cascade. Lateral separation of the instruments in the subcanopy trunkspace has little effect on the shape of the velocity spectra. On the other hand, lateral and vertical velocity coherences between spatially separated sensors are low inside the canopy. These low coherences are due to the Eulerian length scales being of the same order of magnitude as the separation distances of the anemometers. Phase angles between velocity components are about zero for small separation distances. When the two instruments are separated by 9 m and one instrument is positioned in a row while the other is between two rows, vertical velocities are about 180 deg out of phase and the streamwise velocities are about 40 to 60 deg out of phase. These data support the contention that preferred differences occur between within- and between-row wind flow regimes.     

Scale Properties of Anisotropic and Isotropic Turbulence in the Urban Surface Layer

The scale properties of anisotropic and isotropic turbulence in the urban surface layer are investigated. A dimensionless anisotropic tensor is introduced and the turbulent tensor anisotropic coefficient, defined as C, where \(C = 3d_{3}\,+\,1 (d_{3}\) is the minimum eigenvalue of the tensor) is used to characterize the turbulence anisotropy or isotropy. Turbulence is isotropic when \(C \approx 1\), and anisotropic when \(C \ll 1\). Three-dimensional velocity data collected using a sonic anemometer are analyzed to obtain the anisotropic characteristics of atmospheric turbulence in the urban surface layer, and the tensor anisotropic coefficient of turbulent eddies at different spatial scales calculated. The analysis shows that C is strongly dependent on atmospheric stability \(\xi = (z-z_{\mathrm{d}})/L_{{\textit{MO}}}\), where z is the measurement height, \(z_{\mathrm{d}}\) is the displacement height, and \(L_{{\textit{MO}}}\) is the Obukhov length. The turbulence at a specific scale in unstable conditions (i.e., \(\xi < 0\)) is closer to isotropic than that at the same scale under stable conditions. The maximum isotropic scale of turbulence is determined based on the characteristics of the power spectrum in three directions. Turbulence does not behave isotropically when the eddy scale is greater than the maximum isotropic scale, whereas it is horizontally isotropic at relatively large scales. The maximum isotropic scale of turbulence is compared to the outer scale of temperature, which is obtained by fitting the temperature fluctuation spectrum using the von Karman turbulent model. The results show that the outer scale of temperature is greater than the maximum isotropic scale of turbulence.     

Stratified airflow over one or two hills

Methods for calculating, interpolating and idealising air flow in complex terrain are reviewed. Then the general structure of stratified airflow over a single hill of height H and length L₁ is studied in detail and shown to be determined by the upwind velocity profile, the magnitude of a characteristic Froude number and the dimensions of the hill. Let N(L₁) be the buoyancy frequency upwind at a height L₁, and u_* and U₀ be the upwind friction and mean velocity respectively; then the flow is effectively neutral if u^*/NL₁>1. But if u_*/NL₁>1 and u₀/NL₁>1, the stratification is weak enough to affect the upwind turbulence and velocity profile but not the dynamics of the flow over the hill. If U₀/NL₁ <1 but U₀/NH>1 the buoyancy forces are strong enough to affect the mean flow over the hill but not strong enough to prevent it passing over the top. In this regime the flow is very sensitive to the form of the upwind temperature profile. If U₀/NH>1, much of the flow passes round the hill. A similar classification, with different flow patterns, is appropriate for unstably stratified flows. When the wind is weak enough, local slope winds can dominate. Results from the analysis of these different regimes are described and compared with laboratory of field measurements where possible. It is shown how some of these results can be extended to groups of hills.     

Coherent eddies and turbulence in vegetation canopies: The mixing-layer analogy

This paper argues that the active turbulence and coherent motions near the top of a vegetation canopy are patterned on a plane mixing layer, because of instabilities associated with the characteristic strong inflection in the mean velocity profile. Mixing-layer turbulence, formed around the inflectional mean velocity profile which develops between two coflowing streams of different velocities, differs in several ways from turbulence in a surface layer. Through these differences, the mixing-layer analogy provides an explanation for many of the observed distinctive features of canopy turbulence. These include: (a) ratios between components of the Reynolds stress tensor; (b) the ratio K _H/K _M of the eddy diffusivities for heat and momentum; (c) the relative roles of ejections and sweeps; (d) the behaviour of the turbulent energy balance, particularly the major role of turbulent transport; and (e) the behaviour of the turbulent length scales of the active coherent motions (the dominant eddies responsible for vertical transfer near the top of the canopy). It is predicted that these length scales are controlled by the shear length scale % MathType!MTEF!2!1!+-% feaafeart1ev1aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn% hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr% 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq-Jc9% vqaqpepm0xbba9pwe9Q8fs0-yqaqpepae9pg0FirpepeKkFr0xfr-x% fr-xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamitamaaBa% aaleaacaWGtbaabeaakiabg2da9iaadwfacaGGOaGaamiAaiaacMca% caGGVaGabmyvayaafaGaaiikaiaadIgacaGGPaaaaa!3FD0!\[L_S = U(h)/U'(h)\] (where h is canopy height, U(z) is mean velocity as a function of height z, and % MathType!MTEF!2!1!+-% feaafeart1ev1aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn% hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr% 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq-Jc9% vqaqpepm0xbba9pwe9Q8fs0-yqaqpepae9pg0FirpepeKkFr0xfr-x% fr-xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGabmyvayaafa% Gaeyypa0JaaeizaiaadwfacaGGVaGaaeizaiaadQhaaaa!3C32!\[U' = {\rm{d}}U/{\rm{d}}z\]). In particular, the streamwise spacing of the dominant canopy eddies is _x = mL _s, with m = 8.1. These predictions are tested against many sets of field and wind-tunnel data. We propose a picture of canopy turbulence in which eddies associated with inflectional instabilities are modulated by larger-scale, inactive turbulence, which is quasi-horizontal on the scale of the canopy.     

Spectral Short-circuiting and Wake Production within the Canopy Trunk Space of an Alpine Hardwood Forest

Using synchronous multi-level high frequency velocity measurements, the turbulence spectra within the trunk space of an alpine hardwood forest were analysed. The spectral short-circuiting of the energy cascade for each velocity component was well reproduced by a simplified spectral model that retained return-to-isotropy and component-wise work done by turbulence against the drag and wake production. However, the use of an anisotropic drag coefficient was necessary to reproduce these measured component-wise spectra. The degree of anisotropy in the vertical drag was shown to vary with the element Reynolds number. The wake production frequency in the measured spectra was shown to be consistent with the vortex shedding frequency at constant Strouhal number given by f _vs = 0.21ū/d, where d can be related to the stem diameter at breast height (dbh) and ū is the local mean velocity. The energetic scales, determined from the inflection point instability at the canopy–atmosphere interface, appear to persist into the trunk space when , where C _du is the longitudinal drag coefficient, a _cr is the crown-layer leaf area density, h _c is the canopy height, and β is the dimensionless momentum absorption at the canopy top.     

A Theoretical Model for the Study of Convective Turbulence Decay and Comparison with Large-Eddy Simulation Data

The development of a theoretical model fora decaying convective boundary layeris considered. The model relies on thedynamical energy spectrumequation in which the buoyancy andinertial transfer terms are retained,and a closure assumptionmade for both. The parameterization for thebuoyancy term is given providing a factorizationbetween the energy source termand its temporal decay. Regarding the inertialtransfer term a hypothesis ofsuperposition is used to describe theconvective energy source and time variationof velocity correlation separately.The solution of the budget equation for theturbulent kinetic energy spectrum is possible,given the three-dimensional initial energyspectrum. This is doneutilizing a version of the Kristensen et al.(see Boundary-Layer Meteorol. 47, 149–193)model valid for non-isotropic turbulence. During thedecay the locus of the spectralpeak remains at about the sameposition as the heat flux decreases.Comparison of the theoretical modelis performed against large-eddy simulationdata for a decaying convectiveboundary layer.     

Upstream stagnation points in stratified flow past obstacles

An experimental study has been made of stagnation points and flow splitting on the upstream side of obstacles in uniformly stratified flow. A range from small to large values of Nh/U (where N is the buoyancy frequency, h_m is the maximum obstacle height and U is the undisturbed fluid velocity) has been covered, for three obstacle shapes which are, respectively, axisymmetric, and elongated in the across-stream and in the downstream directions. Upstream stagnation for the first two of these models does not occur until Nh_m/U > 1.05, where it occurs at z ≈ h_m/2. On the central line below this point the flow descends and diverges, and we term this ‘flow splitting’. For the third model (elongated in the downstream direction), stagnation upstream first occurs at Nh_m/U ≈ 1.43, at z ≈ 0. Results for this obstacle are not consistent with the ‘Sheppard criterion’, and this upstream flow stagnation is not apparently related to lee wave overturning, in contrast to flow over two-dimensional obstacles.     

Instability in Lagrangian stochastic trajectory models,and a method for its cure

A number of authors have reported the problem of unrealistic velocities (“rogue trajectories”) when computing the paths of particles in a turbulent flow using modern Lagrangian stochastic (LS) models, and have resorted to ad hoc interventions. We suggest that this problem stems from two causes: (1) unstable modes that are intrinsic to the dynamical system constituted by the generalized Langevin equations, and whose actual triggering (expression) is conditional on the fields of the mean velocity and Reynolds stress tensor and is liable to occur in complex, disturbed flows (which, if computational, will also be imperfect and discontinuous); and, (2) the “stiffness” of the generalized Langevin equations, which implies that the simple stochastic generalization of the Euler scheme usually used to integrate these equations is not sufficient to keep round-off errors under control. These two causes are connected, with the first cause (dynamical instability) exacerbating the second (numerical instability); removing the first cause does not necessarily correct the second, and vice versa. To overcome this problem, we introduce a fractional-step integration scheme that splits the velocity increment into contributions that are linear (U _i ) and nonlinear (U _i U _j ) in the Lagrangian velocity fluctuation vector U, the nonlinear contribution being further split into its diagonal and off-diagonal parts. The linear contribution and the diagonal part of the nonlinear contribution to the solution are computed exactly (analytically) over a finite timestep Δt, allowing any dynamical instabilities in the system to be diagnosed and removed, and circumventing the numerical instability that can potentially result in integrating stiff equations using the commonly applied explicit Euler scheme. We contrast results using this and the primitive Euler integration scheme for computed trajectories in a drastically inhomogeneous urban canopy flow.