Appropriate Momentum Provision for Numerical Simulations of Horizontally Homogeneous Urban Canopies Using Periodic Boundary Conditions

Turbulent flow over urban-like roughness has been numerically studied for various purposes, such as the clarification of turbulent characteristics over rough walls, visualization of turbulent structures around block arrays, and evaluation of urban ventilation and pedestrian winds. In such simulations, a portion of the developing boundary layer is extracted as a numerical domain, assuming periodic boundary conditions in the horizontal direction to reproduce laterally homogeneous rough surfaces. However, the conditions required to drive the airflow by an artificial momentum source uniquely determine the turbulent statistics, which are different from those in developing boundary layers. Therefore, this study presents a new approach for driving the airflow over urban-like block arrays. The new method is based on spatially averaged Navier–Stokes equations to prove the necessity of height-dependent momentum provision. The turbulent flow over a cubical-block array is determined using large-eddy simulations driven by four different momentum sources. Regardless of the driving force, the vertical profiles of the streamwise velocity components are identical. In contrast, slight differences in the vertical Reynolds stress, variances in the velocity components, and turbulence kinetic energy are appropriately reproduced in the new approach. In addition, the budget equations of the turbulent statistics reveal that a change in the vertical Reynolds stress affects the energy production and its redistribution into variance components. The proposed method can contribute to the reproduction of a realistic turbulent flow and provide a theoretical understanding of the momentum provision.


Introduction
Airflow affected by urban geometries has a considerable influence on the comfort and safety of pedestrians in urban spaces. Turbulent flow, particularly owing to the strong shear production due to urban canopies, contributes to environmental problems such as the heat island phenomenon, pollutant diffusion, and gusty winds in urban areas because of the momentum and scalar transport between urban areas and the atmosphere. To examine the impact of urban geometries on turbulent flow formation, transport phenomena, wind environment, and resultant impacts on pedestrians, numerous studies on urban turbulent boundary layers have been conducted using simplified roughness, which are sufficiently elongated arrangements of simplified blocks in the horizontal direction, by means of outdoor measurements, wind-tunnel experiments (WTEs), and computational fluid dynamics (CFD) approaches.
Regarding outdoor measurements, for example, Biltoft (2001) conducted outdoor measurements of the airflows around a model site called Mock Urban Setting Test (MUST), which consists of numerous full-scale containers placed to mimic urban-canopy layers, which mimics an infinitely repeated scenario in the horizontal direction. Inagaki and Kanda (2008) also reported a series of outdoor measurements using an outdoor observation site termed the comprehensive outdoor scale model (COSMO). The COSMO consist of 1.5 m cubes more than fifty times the block height in the streamwise direction to experimentally simulate an urban turbulent boundary layer. This site has been extensively used for various purposes, such as analyzing the scaling method for turbulent statistics (Kanda and Moriizumi 2009), turbulent organized structures (Inagaki and Kanda 2008;Takimoto et al. 2011), and investigation of the effect of urban boundary layers on wind-induced cross ventilation (Hirose et al. 2019(Hirose et al. , 2020. By means of WTEs, the simplified urban canopies have been extensively studied. For example, Cheng and Castro (2002) performed velocity and drag measurements over urbanlike surfaces by employing cubes elongated over four hundred times the height to generate an idealized fully developed urban boundary layer. They studied the universality of the normalized velocity component and turbulence statistics. Accordingly, Cheng et al. (2007) employed identical roughness to determine the roughness length and displacement height for surfaces characterized by packing density. Hagishima et al. (2009) reported datasets of the drag coefficients, velocity profiles, roughness heights, and displacement heights by employing a sufficient fetch of more than one hundred times the building height. Hirose et al. (2022) reported the instantaneous flow fields, statistics, and probability densities of two velocity components in the same simplified urban canopy using time-resolved particle-image velocimetry.
In these studies on outdoor measurements and WTEs, horizontally homogeneous airflow fields were appropriately created to minimize the effect of the streamwise development of the turbulent boundary layers by ensuring a sufficient length of the fetch consisting of simple blocks. Owing to these experimental conditions, the influence of urban geometries on turbulent and aerodynamic properties can be discussed.
Two methods are used to reproduce these conditions in CFDs. One is to reproduce the entire wind tunnel (e.g., Nakajima et al. 2018;Michioka et al. 2019;Okaze et al. 2021), or outdoor measurement sites Dejoan et al. 2010;Oh et al. 2011). Naturally, these methods can reproduce the turbulent flow generated by numerous block elements by both shear and wake turbulent production, resulting in the formation of reasonable turbulent flow fields, even in CFDs. In contrast, the method requires a high computational load because of the large numerical domains that reproduce the entire target domain used in the measurements.
The other method is simplistic; it extracts a portion of the measurement sections owing to the repetitive roughness arrangements in the horizontal directions and horizontally imposes periodic boundary conditions (e.g., Kanda et al. 2004;Kanda 2006;Coceal et al. 2006Coceal et al. , 2007aIkegaya et al. 2017;Ishida et al. 2018). This method can reproduce an infinite group of blocks and ensure a horizontal homogeneity. In addition, this method can reduce the cost of the numerical calculations. For example, Kanda et al. (2004Kanda et al. ( , 2006 analyzed turbulent boundary layers over simplified block arrays with a variety of building densities using large-eddy simulation (LES) and discussed the dependence of the turbulent structure on the geometries. Accordingly, Coceal et al. (2006Coceal et al. ( , 2007a performed direct numerical simulations of the airflow around a cubical roughness array and visualized the turbulent structures formed within the urban turbulent boundary layers. Kawaminami et al. (2018) analyzed the airflow around block arrays to investigate the mutual contribution of wind velocity components and scalar concentration at pedestrian levels in urban canopies. As shown in previous studies, CFDs with periodic boundary conditions also allow us to understand the turbulent characteristics, spatial structures of airflow, and their impact on pedestrians.
However, this method has a critical disadvantage. In other words, the method requires an external force to drive the flow, which determines the vertical profiles of the Reynolds stress. Table 1 presents a brief review of previous studies using periodic boundary conditions for vegetation and urban canopy flow. To the best of our knowledge, only two types of driving forces have been employed for CFDs with periodic boundary conditions: imposing a constant pressure gradient in the streamwise direction and imposing a constant shear stress at the upper boundary of the numerical domain (Table 1). The former is classically termed Couette flow, whereas the latter is known as the Hagen-Poiseuille flow. These two well-studied classical flow representations assume the vertical distribution of the momentum stress is either linear or constant with respect to height. Because turbulent momentum flux dominates in turbulent boundary layers, the unique determination of the vertical momentum stresses in these classical flows indicates that the vertical Reynolds stress in these CFDs is uniquely determined by the driving force, as reported by Watanabe (2004) and Ishida et al. (2018).
In contrast, the actual profiles of the vertical Reynolds stress in WTEs do not exhibit constant or linear trends with respect to height, but that of convex and concave curves with an inflection point regardless of the surface type (e.g., Cheng and Castro 2002 for a simplified urban canopy, Harun et al. 2013 for a smooth surface, and Inagaki et al. 2017 for real districts in an urban area reproduced by CFD). This suggests that the momentum supply via the driving force in the numerical calculations significantly differs from the momentum provision mechanism of the boundary layers, although the urban turbulent boundary layers are formed using a long fetch. As a result, differences in the turbulence statistics and turbulence kinetic energy are expected. However, there are no alternatives to driving the airflow over homogeneous urban canopies, except for these two methods, because the differences in the momentum provisions have not been discussed previously. Therefore, the present study aims at three aspects: proposing an appropriate method for the momentum provision to drive the flow over simplified block arrays with periodic boundary conditions by considering the differences between developing and horizontally homogeneous boundary layers based on the momentum budget equations, comparing the conventional methods with the new approach of the driving force in terms of turbulent statistics and momentum budget equations, and clarifying the direct and indirect effects of the driving force on turbulent generation by detailed analyses of the variances and turbulent kineticenergy budget. In Sect. 2, the theoretical derivation for proposing a new momentum provision is explained; in Sect. 3, the details of the numerical simulations are described; in Sect. 4,

Temporally and Spatially Averaged Equations for Vertical-Momentum Fluxes
To consider the momentum budget in numerical domains with periodic boundary conditions and boundary layers formed over rough surfaces, we assumed the boundary layers described Fig. 1 Schematics of a a numerical domain with periodic boundary conditions and b a developing boundary layer. The quantities of H , δ r , and δ are canopy layer, roughness sublayer, and boundary-layer depths, respectively. L z is the numerical domain height for the situation when the periodic boundary conditions are imposed, and u δ is the temporally and spatially averaged freestream streamwise velocity component in Fig. 1, consisting of a canopy layer (CL) with roughness elements, a roughness sublayer (RSL), and an outer layer. The streamwise, spanwise, and vertical directions are denoted as x, y, and z, respectively, and the corresponding velocity components are u, v, and w, respectively. The coordinates and velocity components are also denoted as x i and u i with i 1, 2, and3, respectively. The heights of the CL, RSL, and boundary layer are denoted by H , δ r , and δ, respectively. If horizontally homogeneous roughness elements infinitely continue in x − y directions, H is regarded as a constant value with respect to x and y, whereas δ r and δ are allowed to change in the x-direction, namely, δ r ≡ δ r (x) and δ ≡ δ(x), respectively.
To derive the spatially and temporally averaged governing equations of the velocity components, we introduced a temporal and spatial average of the quantity φ. The temporal average is defined as: where T represents the average duration. When T is sufficiently long (i.e., T → ∞), φ is assumed constant with respect tot. The deviation from the mean is denoted by φ φ − φ. Similarly, the spatial average of φ is expressed as: We considered a control volume (CV) with a sufficiently large horizontal plane L x L y and minute vertical distance z. Thus, the CV is defined as V L x L y z. Additionally, V F for the integral volume in Eq. (2) represents the volume occupied by the fluid within V when CVs below H are considered. Above CL, V V F anywhere. In Eq.
(2), the division by V is called an extrinsic average. Böhm et al. (2013), Kono et al. (2010a; called this averaging method as a type-2 exclusive average. By contrast, Schmid et al. (2019) referred to the averaging from Whitaker (1999) as a superficial average. In contrast, division by V F is called an intrinsic average (Hiraoka 1993;Finnigan 2000;Kono et al. 2010a, b;Böhm et al. 2013 called this method as a type-1 exclusive average, whereas Schmid et al. 2019 refer to this method from Whitaker (1999). When we employ the intrinsic average to derive the governing equations of canopy flow, we need to include the fluid-to-solid ratio (V F /V ) in the terms below H , whereas the extrinsic average does not include this term. Although either method can be applied, the governing equations are interchangeable and yield the same formulation (Appendix 1). The spatial average is defined as the extrinsic average and is applied to temporally averaged quantities, and the deviation from the temporally and spatially averaged values is denoted as φ φ − φ . Using Eqs. (1) and (2), the temporally and horizontally averaged continuity and Navier-Stokes equations for an incompressible fluid can be written as: where p ≡ P/ρ denotes the kinematic pressure, P denotes the pressure, and ρ denotes the fluid density. The term τ i j represents shear stress and is defined as: where the three terms on the r.h.s. represent the Reynolds stress, dispersive stress owing to the spatial nonuniformity of the temporally averaged velocity components, and molecular viscous stress, which is negligible compared with the Reynolds and dispersive stress terms for the turbulent urban boundary layers targeted in this study. Term f i in Eq. (4) represents the total drag acting on the canopy elements due to the viscous stress and pressure differences between the canopy elements in the kinematic unit. The force acting on the canopy elements are defined as positive. The drag force can be formulated by the commutation errors owing to the interchange between the integral and differential (Appendix 1) when the spatial averaging in Eq. (2) was applied to the pressure and viscous-diffusion terms in the Navier-Stokes equations. Therefore, f i can be written as: where S i represents the i-directional component of the surface vector in V F intercepted by the canopy elements. The surface vector S i is positive when we consider the direction from the fluid to solid objects. Note that S i does not include the bottom surface; therefore, the viscous drag on the bottom surface, or τ 13 (0), which is negligibly small for a deep canopy, is not included in f i . Hence, the total drag force represented by u 2 τ , where u τ is the total friction velocity, can be determined as follows As a scaling friction velocity, we employed u 2 * −τ 13 (z p ), where z p represents the height at which τ 13 becomes the peak value. To differentiate between the total and scaling friction velocities, we introduced u τ instead of u * .
Additionally, the pressure can be decomposed into the spatial deviation part p and the constant streamwise gradient γ x as follows: From Eqs. (4) and (8), the gradient of τ i j can be written as: Equation (9) is identical to the temporally and spatially averaged Navier-Stokes equations, but the expression facilitates the determination of the total momentum flux τ i j , particularly for the vertical momentum flux τ 13 . Appendix 2 explains the treatment of the extrinsic average of the pressure gradient term.
When periodic boundary conditions are imposed, horizontal homogeneous conditions can be rigorously applied because ∂ φ ∂ x ∂ φ ∂ y 0 for any temporally and spatially averaged quantity φ , resulting in the no subsidence condition of w 0. This also implies that we can assume L x , L y → ∞ in Eq. (2), yielding a precise value φ ≡ φ (z). Therefore, the equation to express τ 13 can be simplified as: In contrast, when boundary layers develop in the streamwise direction due to the drag force on the surfaces, the spatial averaging operation in Eq. (2) needs to be considered as a spatial filtering operation at the focal position (x, z) with a finite length of L z , but the spatial averaging operation in the y-direction owing to L y → ∞, results in φ ≡ φ (x, z). Therefore, τ 13 can be expressed as: where p ≈ τ 33 + p (z out ) is derived from the order approximation of Eq. (9) with i 3 (Townsend 1976). Here, z out and p (z out ) are the height and pressure above the boundary layer. Since p (z out ) is constant with respect to x, the pressure term in Eq. (11) is represented as ∂τ 33 /∂ x. Equations (10) and (11) are simply the spatially and temporally averaged Navier-Stokes equations for the streamwise momentum, which can be employed to scrutinize the differences in momentum provision for both cases.

Bulk Momentum Balances
To understand how the total momentum is balanced in the boundary layers for the numerical simulations and experiments, we consider the entire balance of the streamwise momentum by integrating Eqs. (10) or (11) between the bottom and top boundaries.

Periodic Boundary Conditions
From Eq. (10), the momentum balance can be written as: When the flow was driven by a constant pressure gradient in the streamwise direction with a free-slip boundary condition at the top of the numerical domain, we have τ 13 (δ) 0. Hence, the total drag represented by u 2 τ can be related to the source term owing to the pressure gradient (e.g., Kanda et al. 2004;Kanda et al. 2006;Coceal et al. 2006Coceal et al. , 2007aXie and Castro 2006;Kono et al. 2010a, b) as: However, when the flow is driven by a moving lid at the top of the numerical domain (Su et al. 1998(Su et al. , 2000Watanabe 2004), γ 0 and τ 13 (L z ) 0, where L z is the domain height, we obtain: Although these momentum balances under periodic boundary conditions are trivial to derive, the purpose of this study is to propose an appropriate method for accelerating the flow over urban-like surfaces under periodic boundary conditions to reproduce the momentum provision in developing boundary layers observed in experiments. To compare the differences between the numerical simulations and experiments, we used these formulations. Equation (13) indicates that the total momentum absorbed on the rough surface, or the total drag force, balances the total momentum given to the numerical domain as the external body force by the pressure gradient (γ δ). In other words, the condition is identical to the pipe flow driven by the streamwise pressure drop. In contrast, Eq. (14) explains that momentum is provided by the domain-top boundary via surface friction when the moving lid is given as the driving force. Because there is no additional source term for the latter condition, the total momentum provided from the top boundary balances the momentum absorbed at the bottom surface as the drag. These two conditions are considered to be extreme cases when the total momentum is only provided by the pressure gradient or only determined by the top boundary. In realistic experimental scenarios, we need to consider the conditions amid these cases, although either condition is commonly employed to reproduce turbulent flow over rough surfaces (e.g., Kanda et al. 2004;Watanabe 2004;Coceal et al. 2006) and compared them with experiments with sufficient fetch (e.g., Cheng and Castro 2002).

Developing Boundary Layer
In the case of developing boundary layers, we cannot precisely assume horizontal homogeneity of the flow fields, as discussed in Sect. 2.1. Therefore, the total momentum balance must be considered in Eq. (11). By integrating Eq. (11) from the bottom surface, z 0, to the boundary layer top, z δ, we obtain the following balance equation: The suffix δ on variables refers to the values at the top of the boundary layer. At z δ, the total momentum flux is zero, τ 13 (δ) 0. In addition, the surface value of u w does not appear due to the no-slip condition at the bottom surface.
According to Eq. (15), the total momentum absorption to the rough bottom surfaces in the developing boundary layers balances with three momentum source terms on the r.h.s. Term I indicates the momentum provision owing to the streamwise pressure gradient given by the acceleration of the freestream velocity u δ (i.e., u δ d u δ /dx γ in the freestream). Term II is the summation of the momentum sink from the top of the boundary layer owing to w δ > 0 and the momentum source owing to the streamwise deceleration of u . Term III represents the source or sink due to the streamwise change in the turbulent strengths of both u and w. According to the order of approximation, the contribution of Term III is negligible compared to that of Term II (Ikegaya 2022).
To scrutinize the physical interpretation of Term II, we simplify the expression by assuming a zero pressure gradient, γ 0. This allows us to assume that u δ constant with respect to x. In addition, integrating Eq. (3) between the surface and the boundary-layer top can express the vertical velocity component w by u .Therefore, S I I is written as: where θ represents the momentum thickness and is defined as (16) demonstrates that S I I is consistent with the momentum provision term explained by the von Kármán's integral law, which is represented by the streamwise change in momentum thickness. von Kármán's integral equation clearly explains how the relationship between macroscopic quantities such as C d u 2 τ / u 2 δ and θ is expressed based on the momentum balances within the boundary layers. Moreover, the equation is concise for intuitively understanding the momentum balances in the boundary layers. By contrast, the present formulation of Eq. (15) with Eq. (16) provides insight into the physical interpretation of the streamwise change in θ . In other words, the momentum thickness increase indicates that the momentum provision from the boundary layers is due to the integrated decay of the streamwise velocity component in the streamwise directions.
Although the streamwise decay of the mean streamwise velocity component becomes small if a sufficient fetch of rough surfaces is employed in the experiments, the formulations in Eqs. (13-15) imply that the momentum provision mechanisms in experimental conditions are different from those in numerical simulations, where the flow is driven by a body force or top boundary under periodic boundary conditions. Conversely, in the experimentally formed boundary layers, it is essential that the streamwise decay of the mean streamwise velocity components is the genuine source of the momentum absorbed at the bottom surface.

Local Momentum Balance
The differences in the total momentum provision mechanisms also affect the vertical momentum fluxes, indicating that the change in the vertical momentum flux profiles is due to the difference in the balanced local momentum. To clarify this aspect, this section explains how momentum fluxes are related to momentum provision.
In either case discussed in Sect. 2.2, the budget equation of the vertical momentum flux can be rewritten as: where s 1 (z) denotes the generalized momentum source term. The term s 1 depends on how momentum is provided to the system. For example, s 1 γ for a pressure-driven flow and s 1 0 for a lid-driven flow. By integrating Eq. (17) between z 0 and an arbitrary height of z, the momentum flux at z, τ 13 (z), is written as: where S 1 and F 1 are defined by the following integrals: From this definition, term S 1 indicates the accumulated momentum source from z 0 to z, and F 1 represents the total drag force between z 0 and z. The total drag force F 1 has a profile only when z ≤ H , and above this height, the term denotes the total drag force, that is, Assuming that the vertical profile of F 1 is identical regardless of the momentum provision method, Eq. (18) indicates that τ 13 is only determined by S 1 . In other words, a momentum provision method automatically determines the vertical profile of τ 13 as a given condition. It should be noted that the impact of S 1 or s 1 on τ 13 becomes marginal near the canopy layer when δ/H is large (namely, deep boundary layers) because the impact of the term z 0 s 1 (x, ζ )dζ is insignificant with respect to F 1 as can be seen in the minor difference of τ 13 in Fig. 2. For example, we consider two ideal cases of the flow driven by a pressure gradient or a moving lid. As for the former condition, τ 13 (H ) −F 1 (H ) + u 2 τ H /δ, while the latter gives τ 13 (H ) −F 1 (H ). Therefore, the values of τ 13 near the canopy top is nearly identical as F 1 (H ) regardless of the momentum source because u 2 τ H /δ becomes sufficiently small (e.g., the effect of the source term is only 10%, 5%, 3% if δ/H is 10, 20, and 30, respectively).
To clarify the relationship between s 1 , f 1 , S 1 , F 1 , and τ 13 , Fig. 2 shows schematics for a lid-driven, pressure-driven, and developing boundary layer. Although the first two cases are very primitive conditions known as the Couette or Hagen-Poiseuille flow, we include them to contrast the difference with a developing boundary layer. Although the canopy profiles are minor aspects of this discussion, the schematics of f 1 are determined by assuming f 1 C d0 u 2 and u u H ex p(a(z/H − 1)) from the simplified canopy profile proposed by Macdonald (2000), where C d0 u 2 τ /u 2 H is the drag coefficient, a 9.6λ f and λ f are the attenuation coefficient and frontal area index of the canopy, respectively, and u H is the wind speed at H . As shown in Fig. 2, the schematic relationship between the source and flux terms, λ f 0.2, was arbitrarily selected. In Fig. 2a, for the lid-driven flow, the momentum is provided only from the top of the domain, and there is no other source within the numerical domain. This implies that s 1 0 and S 1 0, yielding to the formation of a constant flux layer (CFL) with τ 13 u 2 τ . In contrast, in Fig. 2b, the flow is driven by a constant pressure gradient with s 1 γ constant at any height. Consequently, S 1 γ z and τ 13 (z) are proportional to z.
For a developing boundary layer, we are not sure how the local source term s 1 can be expressed owing to the streamwise change in the streamwise velocity component and turbulent intensities. However, numerous previous studies have shown typical vertical profiles of the vertical Reynolds stress, u w (~τ 13 above the RSL for a fully turbulent boundary layer because of the negligible effect of the molecular diffusivity). For example, Raupach (1981) reported outer layer similarity in the vertical profiles of the vertical Reynolds stress, u w , for airflows over several types of surfaces including a smooth wall and rough surfaces based on wind-tunnel experimental data. The vertical profiles show a convex-to-concave curve from lower to higher positions. As for the urban turbulent boundary layers, Inagaki et al. (2017) Fig. 2 Schematics of difference in driving force and source term in the momentum equation, and resultant momentum flux (a1, b1, c1). Schematics of boundary layer (a2, b2, c2), vertical profiles of source and drag terms (a3, b3, c3), vertical profiles of integrated source and drag terms, and resultant vertical momentum flux performed a numerical simulation for the airflow around realistic buildings in a full scale and showed that the vertical Reynolds stress has a convex shape with respect to z, and changes into a concave curve when gradually approaching zero from the top of the boundary layer. Based on these observations, we can qualitatively describe the schematic profile of τ 13 with convex and concave curves from the bottom to the top of the boundary layer (a more detailed determination based on previous experimental data is explained in Sect. 3.2), as shown in Fig. 2c. Consequently, S 1 τ 13 + F 1 and s 1 d S 1 /dz are determined from τ 13 . Because s 1 is determined by the slope of S 1 , we can obtain a bell-shaped curve with a peak value at the inflection point of τ 13 (Fig. 2c). Although these profiles are simple estimations that may depend on the surface conditions as well as the Reynolds numbers, Raupach (1981) showed that such inflection points in u w are universally shown for both smooth and rough surfaces. This means that the momentum source term owing to the mean-velocity decay and turbulence intensity change in the streamwise direction has peak values within the boundary layers, as schematically shown in Fig. 2c.
The schematics of s 1 , S 1 , and τ 13 in Fig. 2 are an example of a qualitative discussion; however, the relationship in Fig. 2 highlights that height-dependent values of the momentum provision are required for reproducing τ 13 observed in developing boundary layers, which cannot be simulated by Couette-or Hagen-Poiseuille-type driving forces. Smaller source values near the bottom surface and top of the boundary layer are plausible because the sources of the developing boundary layer are determined by the terms on the r.h.s. of Eq. (15), indicating that the streamwise change in the statistics becomes milder in such regions compared with the values in the middle of the boundary layer. In addition, the comparison among the three conditions yields the speculation that Couette-and Hagen-Poiseuille-type driving forces create the boundary layers, which are regarded as extreme conditions, when there are zero or constant sources, whereas such conditions cannot be seen in the entire boundary layer formations. When such conditions can be achieved in a part of a boundary layer, the CFL with zero local momentum sources or linear decrease in the momentum flux with a constant local momentum source can be seen.
To summarize this discussion, derivations of the relationships between the source term and resulting momentum-flux profiles enable us to understand how a developing boundary layer can be mimicked using periodic boundary conditions. In other words, the momentum provision method automatically determines the vertical profile of τ 13 as a given condition in a numerical simulation. Therefore, if we provide certain profiles of the momentum source s 1 , which are determined by experimental data or targeted conditions, as an additional term in Navier-Stokes equations, τ 13 profiles can theoretically be reproduced in arbitrary shapes in entire boundary layer. It should be noted that the present method does not deny the reproduction of τ 13 within the roughness sublayer or canopy layer by the CFDs driven by either constant-body force or moving-lid (e.g., Xie et al. 2004Xie et al. , 2006Coceal et al. 2007a, b). This is because the relative impact of the determination of s 1 and S 1 is marginal when the boundary layer is deep or when we consider the lower part of the entire boundary layer.

Numerical Descriptions
LESs were performed for airflows around a simplified urban-block array using OpenFOAM v1912, which is an open source CFD package based on the finite volume method (OpenFOAM 2016). The governing equations are the incompressible filtered continuity and Navier-Stokes equations, with a turbulence model expressed as follows: hereafter variable u i represents the grid-scale velocity components of the LESs. The subgrid stress τ SG S i j ν SG S ∂u i /∂ x j is modelled using the standard Smagorinsky model with , and x i is the grid size in the i-direction. We adopted the van Driest damping function to reduce ν SG S near the wall. The value p indicates the spatial deviation of the kinematic pressure and s i is the driving source term. The temporal and spatial derivative terms were discretized using second-order backward scheme and the second-order linear interpolation scheme, respectively. The time step t in the temporal development was set to 10 −3 s. The unsteady SIMPLE (semi-implicit method for pressure linked equations, Versteeg and Malalasekera 1995) method was employed to couple velocity and pressure. Figure 3 shows a schematic of the calculation domain. The calculation domain was set as 4H × 4H × 4H in the x, y, and z directions, respectively. The building height (H ) was 0.1 m. Four cubes (H × H × H in size) were arranged in a staggered-array layout with a packing density of 25%. A uniform grid resolution of H /20 was employed in each direction based on a previous sensitivity study of grid resolutions to reproduce both turbulent statistics around and above the simplified block array (Coceal et al. 2006). Periodic boundary conditions were applied in both the streamwise and spanwise directions to reproduce an infinitely repeated and simplified urban array with no-slip boundary conditions imposed on the floor of the domain and block walls.
Regarding the numerical settings in the present simulations, we mention two aspects: grid resolution and domain size. The effect of grid resolution on the reproduction of turbulent statistics was evaluated by Xie and Castro (2006) between H /64 and H /8, and they concluded that H /16 to H /64 shows reasonable agreement with each other and experimental data. As for the domain size, although the periodic boundary conditions provide repeated geometrical conditions in both lateral directions, Coceal et al. (2007a, b) demonstrated that larger-scale turbulence is not fully generated when the numerical domain is small, resulting in an insufficient reduction of the two-point correlation coefficients at the boundaries. In contrast, they also showed that the vertical profiles of the mean wind and turbulent statistics were not affected by the domain size and were comparable with the results of the WTEs  (Cheng and Castro 2002). Because our purpose is to propose an appropriate method of driving force to achieve vertical momentum flux, we employed the confined but tested domain size in Coceal et al. (2006).

Simulation Cases and Momentum Source
To confirm the discussion in Sect. 2 and the reproduction of the given momentum flux profiles, four simulation cases with different driving forces were considered, as listed in Table 2. Cases 1 and 2 are the reference cases of the solo lid-driven or pressure-driven flow. Case 1 imposed a constant streamwise wind speed on the upper surface as the boundary condition, and Case 2 imposes a constant streamwise pressure gradient via an external force term, such as s 1 −γ .
The other two cases were designed to reproduce the vertical momentum fluxes given in a referring CFD or experiment in which the developing boundary layers were studied. Inagaki et al. (2017) used lattice-Boltzmann LES for the airflow around a realistic building array of 5 × 20 km in Tokyo, Japan. They showed how the characteristic quantities, such as the boundary layer, displacement, momentum thicknesses, vertical profiles of wind speed, turbulent statistics, and vertical momentum fluxes, undergo changes in the streamwise direction. We selected the vertical momentum flux averaged in the downmost spanwise and 16-19 km streamwise areas as a rough-wall example with a realistic scale. The boundary layer height Reynolds number (Re δ u δ δ/ν), is approximately 3.0 × 10 8 . In another case, Scaling friction velocity u 2 * − uw | peak − uw | peak − uw | peak − uw | peak L z 4H : numerical domain height. − uw | peak indicates the peak values of − uw in the vertical direction Harun et al. (2013) examined the mean velocity profile and turbulent statistics under three pressure gradient conditions for a smooth surface, where Re δ 8.3 × 10 4 . Vertical momentum flux was selected as another example. Although both showed the vertical Reynolds stress (but not the total momentum flux), the total momentum flux was mostly identical at these Reynolds numbers because of the negligible contribution of the molecular diffusion flux. We selected these datasets because the vertical profiles of the Reynolds stress in these studies have different trends to the height owing to the difference in the Reynolds number and surface roughness conditions. (i.e., Harun et al. 2013 employed a smooth surface in a wind-tunnel scale, while Inagaki et al. 2017 performed a full-scale CFD simulation for an urban canopy surface). Therefore, we think these two datasets are suitable to highlight the applicability of the present method.
As can be seen from the various results of the vertical Reynolds stress, their shapes change from convex-to-concave between the bottom surface and the top of the boundary layer. As one of the empirical functions to express the shape, we selected the Gaussian distribution and error function pair for the source term s 1 , and integrated source S 1 as: here μ, σ , and L z are the fitting parameters representing the peak height of the momentum source, the standard deviation, and the domain height, respectively. E is the integrated parameter such that the total friction velocity u τ satisfies the following relation: Figure 4 shows a comparison of Eqs. (23) and (24) with reference data. In the source term of Fig. 4a, c, the experimental data was determined by taking the gradient of the vertical Reynolds stress, namely d −u w /dz from the reference data using the central difference scheme. The vertical total momentum fluxes and source terms of Cases 1 and 2 are also drawn for comparison. As can be seen in the figures, the regressed lines in Fig. 4b, d agree with the reference data, except for the lowest regions where dispersive, canopy drag, or viscous drag become significant. Accordingly, the source term in Eq. (23) (the regressed lines) in Fig. 4a, c mostly agreed with the reference data (red lines). Near the wall regions, the values of S 1 , corresponding with the integrated source term, are required to approach u 2 τ , whereas −u w approaches zero. Hence, the extrapolation of S 1 and s 1 based on Eq. (23) is plausible, rather than simply determining the value of d −u w /dz. Near the boundary layer top, when we have available data on the vertical Reynolds stress, we can determine the momentum source in the entire boundary layer to exactly reproduce the shape. However, we need to use the gradient of the momentum fluxes, resulting in possible numerical errors to determine the gradient from the small values of −u w ; thus, we employed Eq. (23) as an empirical equation that represents the source term under these two conditions. Based on the above discussion, these two cases are referred to as Cases 3 and 4, respectively, and mimic the conditions for developing boundary layers. The magnitudes of Cases 3 and

Post Processing
The time duration for the data sampling used in statistical calculations was 275 s for Case 1 and 100 s for Cases 2 to 4, respectively, after an initial calculation duration of 200 s to eliminate the effect of the initial condition. These times satisfy the criteria values of 200T and 400T , where the reference time T is defined by H /u τ , for the initial and sampling calculation periods proposed by Coceal et al. (2006), respectively. In Case 1, the wind speeds near the CL were low because the flow was driven by the top smooth surface as a moving lid. Therefore, a long duration is required for statistical convergence. Turbulent statistics were calculated at 40 × 40 × 40 points in each direction using the time-series data sampled at 1000 Hz.
The calculation results were normalized by the friction velocity u τ or the scaling friction velocity u * , depending on the purpose of the scaling. The scaling friction velocity u * is defined by the peak value of (− uw ) 0.5 . In Case 1, because the CFL is formed, u τ u * , whereas the other three cases yield u * < u τ .
The validations of the mean streamwise velocity components are shown in Appendix 3 by comparing the present data with those obtained by previous LES (Ikegaya et al. 2016) and WTEs (Cheng and Castro 2002), and they show good agreement with each other. Because our simulation conditions follow the many previous sensitivity studies of the grid resolution and numerical domain size, we focused on the post-processing data analyses of the turbulent flow fields in this study.

Reproduction of Vertical Reynolds Stresses
Before starting the comparison of each case, this section justifies the theoretical determination of the momentum fluxes using the given source term, as in Eqs. (18)-(20). Figure 5 presents a comparison of the horizontally averaged Reynolds stresses of Cases 3 and 4, reference data (Inagaki et al. 2017;Harun et al. 2013), and theoretical lines determined in Eq. (19). These values were normalized by the total friction velocity u τ to confirm the ratio of the terms to the total source term given to the numerical domain.
In Fig. 5a, b, both the vertical Reynolds stresses, and integrated source terms have similar tendencies above the CL, showing a mild upward and downward convex curve. In addition, the slight differences in the vertical Reynolds stress of Cases 3 and 4 were reproduced by providing different source terms. In addition, either of the profiles in Cases 3 and 4 is considerably different from those in Cases 1 and 2 (Fig. 2), which are theoretically determined as a linearly decreasing line or constant value. According to these results, the given source term, which mimics the momentum source owing to the velocity decay in the developing boundary layer, can precisely reproduce the target profiles of a vertical Reynolds stress.
In contrast, in the CL, the integrated source terms and reference data agreed with the results of the parameter fitting. However, the calculation results were different for these two lines in each case. This is owing to the effect of the drag force acting on the blocks (Eq. 17) and dispersive stress, which are constituents of the total momentum in the CL. In the present methods, which provide arbitrary momentum sources, it is possible to reproduce the total momentum fluxes in agreement with the target shape. However, the vertical Reynolds stress within the CLs was determined because of the force on and airflow around the canopy. Therefore, we need to validate the vertical Reynolds stress or drag force within the CLs based on the experimental data conducted under the same geometrical conditions. The deviation in the present CFD and reference data in the CL is due to the difference in CL heights; therefore, we cannot directly compare the values within the CL.

Mean Wind Speed and Turbulent Statistics
A major concern in applying the periodic boundary conditions is whether previously used methods such as Cases 1 and 2 are different from those determined by Case 3 or 4, in which at least the vertical Reynolds stress agrees very well with the reference data. Therefore, we investigated the effect of the momentum source on the reproduced mean wind speed and turbulence statistics. Figure 6 presents a comparison of the profiles of the horizontally and temporally averaged streamwise velocity components u for the four cases. The velocity component was normalized by U 2H , which is defined as the horizontally and temporally averaged streamwise velocity at z/H 2.0. As shown in the figure, there are no substantial differences in the profiles among Cases 1 to 4 in the range of 0 < z/H < 3.0. Above z/H 3.0, the profile of Case 1 showed a clear difference because the flow was driven by the smooth wall at the top, which created a smooth turbulent boundary layer at the top of the numerical domain. The agreement between Case 1 and 2 was reported in previous studies (e.g., Watanabe 2004) if appropriate normalization was applied. In addition, there were no apparent differences in the profiles of Cases 3 and 4, implying that the driving force of the flow was less important in determining the mean velocity fields. This fact justifies numerous previous studies on urban boundary layers conducted under periodic boundary conditions; constant-pressure Fig. 6 Vertical profiles of horizontally and temporally averaged streamwise velocity component u for Cases 1 to 4. U 2H is the horizontally and temporally averaged wind speed at z 2H driven, lid-driven, or flow driven by vertically changing sources does not cause substantial differences.
The difference in the momentum source must be shown in other statistics, as we expected a difference in vertical momentum flux. Figure 7a, b show the vertical profiles of the horizontally averaged Reynolds stress − u w and dispersive stress, − u w normalized by u 2 τ . For LESs, the vertical momentum flux in Eq. (5) is added to the subgrid scale (SGS) viscous term ν SG S ∂u/∂z . The viscous term ν∂ u /∂z and the SGS viscous term are negligibly small because of the sufficiently large Reynolds number and fine grid resolution. Therefore, we only showed − u w and − u w The height where the Reynolds stress becomes peak and the dispersive stress is nearly equal to zero was defined as the RSL depth.
In Fig. 7a, the profiles from Cases 1 to 4 exhibit markedly different tendencies above the RSL. In Case 1, the Reynolds stress is approximately constant at a height between 1.2H and 3.5H . In contrast, the profiles of Cases 2 to 4 show a decreasing trend with increasing height. The trend is linear in Case 2, whereas the slopes of Cases 3 and 4 with respect to height are slightly curved. Owing to the horizontal homogeneity, the dispersive stress in Fig. 7b is slightly above the RSL. These differences in the vertical profiles of u w is simply owing to −∂ u w /∂z s 1 (z). Hence, the driving force represented by s 1 (z) directly affects the vertical gradient of the Reynolds stress, as shown in Fig. 7a.
In addition to the differences in the vertical profiles of the Reynolds stress, Fig. 7 shows interesting characteristics of the scaling of the Reynolds and dispersive stresses. In the RSL, the profiles of Reynolds stress from Cases 1 to 4 have different values, as shown in Fig. 7a below z 1.2H , whereas the profiles of dispersive stress agree well regardless of the difference in the momentum source, as shown in Fig. 7b. In contrast, the Reynolds stress below 1.2H is consistent when u 2 * is employed as the scaling velocity instead of u 2 τ (Fig. 7c). These scaling characteristics indicate that the dispersive stress is universally scaled by the total momentum source denoted by u τ , which determines the mean velocity components because it represents the spatial heterogeneity of the mean wind velocity component distributions. In contrast, the Reynolds stress within the RSL and CL is influenced by the peak values of the total momentum flux τ 13 (z p ) because u * is the representative variable of the temporal perturbation of the velocity components. This implies that defining a scaling friction velocity by the peak values of the Reynolds stress is a plausible method for WTEs, CFDs, and field measurements because we can consider the differences in the driving force of the flow by employing u * . However, the total momentum source is not always scaled with this velocity, particularly when the flows are driven by different momentum sources as depicted in Fig. 7a.
The vertical profiles of the horizontally averaged turbulent kinetic energy (TKE), e t 0.5 u 2 i where e t ≡ 0.5u i 2 , and the variance of the streamwise, spanwise, and vertical velocity components, u 2 , v 2 , and w 2 are presented in Fig. 8. The TKE and variances were normalized by scaling the friction velocityu 2 * . In Fig. 8a-d, all the profiles agree well with each other below the RSL and show a concave shape with respect to height. This means that the ratio of the TKE or variances to u 2 * is identical regardless of the case, implying that the turbulent flow below the RSL is not affected by the driving force.
In contrast, the profiles of the TKE and variances above the RSL were different among cases. In Case 1, the profiles were approximately constant in the vertical direction below 3H , whereas those of the other cases decreased with height.
Although there were slight differences in the profile shapes among Cases 2 to 4, all of them showed remarkably similar tendencies against the height. In other words, the values of the TKE and variances of u and v peaks at the top of the RSL show an upward convex below 2H , and become concave above this height. The variances of w show peak values slightly above the top of the RSL and decrease linearly with height.
As can be seen in the differences in the momentum flux τ 13 directly influenced by the driving force, it is thought that the differences in the TKE and variances occur because of the indirect influence of the driving force on TKE production.

Budget of the Turbulent Kinetic Energy
To consider the reason for these differences in TKE among the cases, the budget equation of the spatially averaged TKE is discussed in this section. The budget equation for e t is written as: represent the TKE production, turbulent transport, pressure transport, molecular diffusion, and dissipation, respectively, determined by the resolved scale velocity components. The contributions of the SGS turbulence are represented by the SGS transport, Figure 9 shows the vertical profiles of the horizontally averaged terms in Eq. (26). These terms are normalized by u 3 * /H . The residual part includes molecular dissipation by the resolved velocity components, SGS transport, and SGS production terms. The SGS transport term was not quantified in our simulations because its contribution to the residual is thought to be small. Therefore, we can assume that the residual is an approximate summation of the other two terms. Because we employed an LES, the total dissipation rate of the TKE could not be precisely determined, and the physical meaning of the residual must be considered as the total transport of the turbulent energy into the smaller than the SGS turbulence. Therefore, The two transport terms, T k t and T p t , exhibit similar tendencies in all cases. In other words, the transport terms are convex with positive values in the CL, reach a negative peak at approximately 1H , and tends to zero from above the RSL height. These are well-known characteristics of turbulent transport in the CL, where energy produced near the canopy top owing to the strong velocity shear at 1H is transported to the lower level of the canopy. As can be seen in these profiles, the contributions of the driving force to these energy transports are marginal, probably because the spatial energy transport is characterized by cubes but not by the driving force applied to the flow.
Similarly, the differences in the diffusion term D t v between the cases were also negligible. In addition, the ratio of the term to the residual and production terms was significantly small. A notable effect of the driving force can be confirmed in both the production and residual terms. The production term of Case 1 (Fig. 9a) shows a different tendency with respect to the height from those of Cases 2 to 4. Moreover, there were slight differences in the shapes of the production terms among Cases 2 to 4. The same tendency was observed for the residual terms. The residual terms are calculated from the other terms, assuming the budget of the TKE. In addition, both the transport and diffusion terms were hardly affected by the driving force. Therefore, we can conclude that the major effect of the driving force appears in the production term, resulting in a change in the residual (or energy transport into the SGS energy) terms.
To scrutinize how the driving force affects the vertical profile in the production term, Fig. 10a compares the vertical profiles of the production terms normalized by u 3 * /H for the four cases. In the CL, these profiles show the same trends among the cases, with downward convex shapes and a peak near the top of the RSL with a distinct increase near 1H . Although the general trends were remarkably similar among the four cases, only Case 1 showed a milder peak than the other three cases. The production terms decrease with height in the range of 1.2 < z/H < 1.6 with a similar slope. However, above z/H 1.6, the profile in Case 1 departs from the other three cases and changes from the decrease to increase as the height approaches the top of the numerical domain. In contrast, the production terms of Cases 2 to 4 continuously decreased with height and approached zero. Although these general trends in the production terms are similar, there are slight differences in the values of the production terms among the three cases.
To clarify the effect of the roughness elements and velocity shear, we considered the breakdown of P t . The term is decomposed into two terms using horizontal averaging. Because periodic boundary conditions are assumed in the present simulations, the production terms are written as: where P s − u w ∂ u ∂z and P w − u i u j ∂u i ∂ x j represent shear and wake production, respectively (Finnigan 2000). The shear production term refers to the amount of kinetic Fig. 10 Vertical profiles of horizontally averaged a total TKE production term P t and b the breakdowns of the shear and wake production terms, P s and P w , respectively. The horizontal dot dashed line indicates the roughness sublayer height energy transported from the mean kinetic energy to the TKE, defined by the product of the horizontally averaged Reynolds stress and vertical gradient of the horizontally and temporally averaged streamwise velocity. On the other hand, wake production indicates the amount of kinetic energy from the dispersive kinetic energy (DKE) to the TKE because of the spatial heterogeneity of the velocity component distributions. Finnigan (2000) explained that the wake production for vegetation canopy layers is due to the spatial distribution of the flow fields with the scale of roughness elements. Similarly, Coceal et al. (2006) showed that the spatial heterogeneity of the velocity components in the CL of cubical arrays requires the consideration of wake production. By contrast, the DKE and wake production terms are not necessarily considered for a horizontal homogeneous flow. Figure 10b shows the vertical profiles of the shear and wake production terms normalized by u 3 * /H . In the RSL, the shear and wake production terms exhibited similar tendencies among the cases. Below 0.8H , the values of P w were larger than those of P S . This is probably because the spatial deviations of the flow fields due to the roughness elements are significant at the lower level of the CL. In contrast, the Reynolds stress, which consists of the temporal deviations of the velocity components, becomes weak owing to the large form drag acting on the roughness elements. In the range of 0.8 < z/H < 1.2, the wake production terms are negative, and the shear production terms steeply increase with the height. The steep increase in P s was due to the strong shear at 1H, whereas the negative values of P w implied that the energy produced at the shear layer was transported to the lower level of the CL via the wake production term. As can be seen, the differences in the energy production mechanisms in the RSL are identical regardless of how the flow is driven.
The effect of the driving force became apparent above the RSL. As the spatial deviations in the Reynolds stress and mean wind velocity in the streamwise velocity component were negligible owing to the horizontal homogeneous flow above the RSL, the values of P w became zero, as shown in Fig. 10b. Therefore, notable differences in the energy production terms can be observed in the vertical profiles of P s in Fig. 10b, resulting in the change in P t shown in Fig. 10a.
Owing to the decomposition of P t into P w and P s , we can clarify the effect of the driving force via the Reynolds stress, because P s is the product of u i u j and the velocity gradient. In Case 1, the Reynolds stress is approximately constant, and the vertical gradient of the streamwise velocity component increases dramatically near the top of the numerical domain owing to the boundary condition at the upper surface. In Cases 2 to 4, the values of u i u j monotonically decrease with height and approach zero near the top of the domain, although the values differ slightly from each other (Fig. 6). In contrast, the vertical profiles of the streamwise mean velocity component were identical (Fig. 5), indicating that similar vertical gradients could be formed for Cases 2 to 4. According to these differences, the shear production terms, represented by the product of these two terms, vary sensitively with differences in the driving force. As a result, the vertical profiles of the TKE and variances differ owing to TKE production mechanisms.

Budgets of Variance
Another concern is how these differences in TKE production owing to the driving force affects the variances of each velocity component. Hence, we discuss the budget of the variances u 2 α of each velocity component u α based on the TKE. Here, α is defined as an index not taking the summation. The budget equations of u 2 α are written as: x j represent the total production, redistribution, turbulent transport, pressure transport, molecular diffusion, and dissipation terms, respectively, as determined by the resolved scale velocity components. The contributions of the SGS turbulence are represented by the SGS transport T SG S Figures 11,12 and 13 shows the vertical profiles of the horizontally averaged terms in the variance budget equations of the streamwise, spanwise, and vertical velocity components for the four cases. These terms were normalized by u 3 * /H . Similar to the TKE budget equation, the residual part can be interpreted as energy transport into the SGS turbulent energy. The term T p u α appears only in the vertical velocity component because of the horizontal average (i.e., ∂ u p /∂ x 0 and ∂ v p /∂ y 0). For all the velocity components, the vertical profiles of the transport terms, T k u α and T p u α , were barely affected by the difference in the driving force. In addition, the production terms for P v and P w were comparable among the four cases. Furthermore, the diffusion term D V u α are negligible. Therefore, in the variances of each velocity component, the effect of the driving force distinctly appears in the production term of P u , the redistribution terms of R u , R v , and R w , and the residuals. This implies that the differences in the momentum source do not change the spatial distributions of the turbulent quantities, whereas it changes the vertical profiles in the Reynolds stress, yielding a change in other terms of the production and redistribution and the resultant energy transfer rate into the SGS turbulence. In Sect. 4.5, we investigate how the driving forces change the production and redistribution terms.

Energy Production and Transfer Among Variances
Similar to the TKE production term, the production terms of the variances (P u α ) are decomposed into two terms, by taking the horizontal average as follows: Here P s represent the shear and wake production of the variances, respectively. Because periodic boundary conditions are applied in the present simulations, the shear production terms of the spanwise and vertical components, P S v and P s w , are zero owing to v w 0 and ∂ φ /∂ x ∂ φ /∂ y 0. The other production terms are expressed as P s u −2 u w ∂ u ∂z , indicting that shear production only works on the streamwise component, whereas wake production remains below the RSL for each velocity component. Figure 14a shows the breakdown of the production terms, namely, P s u , P w u , P w v , and P w w for all cases. By definition, P s u 2P s ; therefore, the vertical tendency of P s u is the same as that of the TKE in Fig. 10b. In the RSL, similar to the breakdown of P t of the TKE, the shear and wake production terms are nearly identical in all cases. Below 0.8H , the values of P s u and P w u i are marginal. However, P s u and P w u become positively and negatively larger between 0.8H to 1.2H , respectively. This tendency is consistent with that in P s and P w for TKE (Fig. 10).
In contrast, above the RSL, the values of P s u differed among the cases because the driving force was the same as those in the TKE, whereas P w u α was negligible because of the horizontal homogeneity of the flow distributions. The reason for the clear differences in P s u is that the driving force directly affects the vertical profiles in u w (Fig. 7) as discussed in terms of P s . To summarize, the effect of the driving force on the variance production terms can only be observed in the production term of the streamwise variance.
To investigate the reason for the difference in the variances of the spanwise and vertical velocity components above the RSL, Fig. 14b shows the vertical profiles of P s u and the redistribution terms of each directional component, R u , R v , and R w (Eq. 28). For both variances of the spanwise and vertical components, the redistribution terms work as the dominant production terms because there is little contribution of wake production to variance production (Fig. 14a), even within the RSL. Above the RSL, the redistribution terms are the only sources because they do not contribute to wake production. Interestingly, the effect of the driving force can be observed in R u , R v , and R w , above the RSL similar to those in P s u . To discuss the relationship between the balance between the shear production and redistribution terms, Fig. 14c shows the ratio of the redistribution of each directional component to the shear production of streamwise variance.
In the RSL, the ratios of the redistribution terms as well as the cases vary with height. Because the vertical Reynolds stress can be scaled by u * but not by u τ (Fig. 7), it is expected that the redistribution mechanisms within the RSL may depend on both local shear production directly affected by the driving force and wake production due to the roughness blocks.
In the range of 1.2 < z/H < 3.0 above the RSL, the ratios of R u , R v , and R w to P s u are approximately constant with respect to height. In other words, the rate of the streamwise component is approximately −2/3, and those of the spanwise and vertical components are 1/3, indicating that the summation of the redistribution is zero. These ratios explain how the TKE is produced, distributed into each velocity component, and eventually dissipated. Above the RSL, the production of streamwise variance was the only source for the TKE because of the horizontal homogeneity. The values of R u /P s u ≈ −2/3 indicates that two-thirds of P s u is transferred into the spanwise and vertical components via R v and R w , and one-third of P s u is dissipated into the internal energy by the dissipation term in the budget equation of the streamwise variances. Because there is no net production in the variances of the spanwise and vertical components, terms R v and R w are the main sources of these quantities. Accordingly, they dissipate in each component. Therefore, we can conclude that the amount of energy loss due to dissipation in each direction is nearly identical because of the isotropic characteristics of small-scale turbulence. We can emphasize that these ratios of redistribution and energy dissipation are the same regardless of the driving force; the change in the driving force only changes the Reynolds stress profiles. However, it indirectly changes the production term in the streamwise variance via the Reynolds stress and the redistribution terms in the spanwise and vertical variances via the production terms, and eventually alters the vertical profiles of the TKE and its breakdown of the variances.

Conclusion
This study proposed a new approach to drive the airflow over urban-like block arrays in CFD by applying periodic boundary conditions to reproduce appropriate vertical Reynolds stresses and TKE based on theoretical approaches that examine the momentum budget equations. The turbulent flow fields over a cubical-block array in a staggered layout with a 25% packing density were simulated using LESs driven by four different momentum provisions. The two cases were determined by the driving force regressed with the previous experiment, along with numerical simulations to specify the upward and downward convex curves of vertical Reynolds stress. The others were based on the conventional two methods, that is, driving the airflow by the constant streamwise pressure gradient and by the moving lid at the top of the domain.
The new method was theoretically derived based on spatially averaged Navier-Stokes equations of airflow over roughness elements. By scrutinizing the differences in the momentum provision in a case with periodic boundary conditions and the situation in which a boundary layer gradually develops in the streamwise direction (e.g., airflow over block arrays by WTEs), we showed that the momentum source in developing boundary layers has heightdependent values. In contrast, neither of the previous methods reproduced these profile types because of the constant momentum provision, resulting in theoretically unique profiles of the vertical Reynolds stress. To compensate for the differences in the momentum provision mechanics, we proposed a method to drive the flow by providing the height-dependent pressure gradient determined by regression using the reference data.
The effect on momentum provision was examined by conducting LESs. Regardless of the differences in the driving forces, the vertical profiles of the streamwise velocity components were identical. The results clearly show that our method can appropriately reproduce slight differences in the vertical Reynolds stress observed in the reference data. In addition to the vertical Reynolds stress, the analyses showed that the driving force affected the variances of the velocity components and TKE.
To clarify how the driving force caused a difference in variances and TKE, the terms of the budget equations were examined. Because of the change in the vertical Reynolds stress, the shear production term was directly affected by the driving force, resulting in a change in the residual term, mostly attributed to energy dissipation into subgrid-scale turbulence. As a result, the vertical profiles of the TKE change owing to differences in the driving force.
Interestingly, the differences in the production term indirectly affected the variance of each velocity component. According to the variance budget equations, the redistribution term in the streamwise component was directly affected by the production term. Because the streamwise redistribution term is the main source for both the spanwise and vertical variances, the difference in this term indirectly affects the turbulent production of the spanwise and vertical components. Despite the differences in the driving force, the ratios between the streamwise production and redistribution terms of each velocity component were nearly identical above the roughness sublayer. This result indicates that the changes in the vertical Reynolds stress owing to the driving force can affect all variances because the turbulent redistribution mechanisms do not change with respect to momentum provision. Therefore, reproduction of the vertical Reynolds stress is essential for reproducing the TKE and variances comparable to those in WTEs.
Although we highlighted the difference in the vertical Reynolds stress, TKE, and variances owing to the momentum source term, it should be noted that this study does not deny previous numerous studies for urban-like canopies employing periodic boundary conditions with a constant pressure gradient as they showed very good agreements of the statistics near the canopy height with experimental data (e.g., Kanda 2006;Coceal et al. 2007a, b;Xie et al. 2004Xie et al. , 2006Kono et al. 2010a, b). The good agreements among these cases are attributed to the fact that the impact of the source term on the turbulent statistics are marginal near the canopy height when we consider a deep boundary layer with large δ/H , which can easily be seen in the full-scale atmospheric boundary layer. In contrast, the boundary layers tend to be shallow for WTEs due to the limitation of facilities, an appropriate source term in the CFDs with periodic boundary conditions needs to be considered.
Although the present method is targeted to boundary layers over urban-like roughness, it can be applied to other simulations where the profiles of the Reynolds stress are specified as a given condition under periodic boundary conditions. Periodic boundary conditions are not a method that perfectly reproduces the conditions for large-scale turbulence, but it is plausible as a simple method to mimic infinitely continuing boundary layers, as employed in numerous previous studies. The proposed method can contribute to such applications by reproducing realistic turbulent flows.
The second term on the r.h.s. of the first line, owing to commutation errors, is zero because u i u j 0 on S k . Similarly, the drag term in Eq. (4) can be derived by applying Eq. (30) to the pressure and viscous terms as: The last three terms on the r.h.s. in the second line appear because of commutation errors, but the last term is zero because u i on S j is zero. The remaining two terms, based on the errors, indicate the pressure-form and surface viscous drag on the canopy elements: These derivations provide the governing equations in Eqs. (3) and (4) based on the extrinsic average. In the extrinsic average, only the boundaries of the obstacles in V must be considered to generate additional terms in the averaged equations. Because the CV of V does not change with respect to x, the interchange between the average and differential is simple.
In contrast to the extrinsic average, the intrinsic average requires consideration of the change in the CV. The intrinsic average of φ, denoted as φ I , is defined as: Owing to the canopy elements, V F ≡ V F (x i ) below z < H . The general form of Eq. (34), taking a convolution with a filter function G(x i ), is written as: here G ≡ G(x i , V F (x i )) 1/V F (x i ) for consistency with the intrinsic averaging filter. The commutation errors are derived by considering ∂ φ I /∂ x i as: The first term on the r.h.s. of the second line arises because of the chain rule of the derivative, and the last term is the commutation error owing to boundaries within V F being similar to the extrinsic average. The clear difference between the extrinsic and intrinsic averages is the additional second term on the r.h.s. of the third line. When the intrinsic average defined by Eq. (34) is applied, the differential in the integrand of the second term can be written as: where m(x i ) V F (x i )/V denotes the fluid-to-volume ratio. The term in Eq. (37) is not a function of x i ; the second term on the r.h.s. of Eq. (37) is: Therefore, the interchange of the average and differential for the intrinsic average is given by the following formula (e.g., Hiraoka 1993): By applying the intrinsic average defined in Eq. (39), the averaged continuity and Navier-Stokes equations can be expressed as: The set of equations is consistent with those derived by Hiraoka (1993) and Kone et al. In Cases 3 and 4, because the vertical profiles of the vertical Reynolds stress are given by the error function in this study, the derivation of the scaling friction velocity cannot be determined algebraically because Eq. (44) is written as: According to Fig. 7a, the Reynolds stresses of Cases 1 and 2 change from the largest to the smallest values, whereas those of Cases 3 and 4 lie between these two cases. In the range of 1.2 < z/H < 3.0. Hence, we consider that the effective friction velocity (u E f f ective * ), in Cases 3 and 4 satisfies: Figure 15 shows a comparison of the vertical profiles of the temporally averaged streamwise velocity component at four different horizontal locations from the present LES and Fig. 15 Comparison of vertical profiles of temporally averaged streamwise velocity component u at four horizontal positions. The black lines represent the present LES results, the red line is the previous LES in the same block arrays as Ikegaya et al. (2016), and the red open circle indicates the WTE results of Cheng and Castro (2002). The values are normalized by the effective friction velocity u E f f ective * reference data from the previous LES (Ikegaya et al. 2016) and WTEs (Cheng and Castro 2002) over the same block array. The velocity component was normalized using the effective friction velocity u E f f ective * . As discussed, the data of Cases 3 and 4 are shown with the possible ranges when u E f f ective * u τ and u τ √ 1 − d/δ. As shown in Fig. 15, the present results for Cases 1 and 2 show good agreement with the reference data. In addition, the results of Cases 3 and 4 normalized by u τ were closer to those of Cases 1 and 2 and the reference data. This is plausible because the vertical slopes of the vertical Reynolds stress in Cases 3 and 4 are milder than those in Case 2. Therefore, we can conclude that the present simulations reproduced the mean velocity fields over the block array.
However, it should be noted that the present results below the RSL, particularly in Fig. 15d, are smaller than those of the WTEs, whereas they agree well with previous LES results (Ikegaya et al. 2016). According to Coceal et al. (2006), the mixing length l m linearly increases with the height above the RSL, where the logarithmic profile is applicable based on the counter-gradient theory in Eq. (44). Therefore, within the RSL and CL, the effective friction velocity is different from u E f f ective * derived herein.