PIV-based pressure estimation in the canopy of urban-like roughness

In-plane velocity measurements from PIV are used to estimate the pressure field above and within the canopy of two staggered arrays of cuboids, with distinct height distributions, via 2D-RANS and 2D-TH. The viability of this approach is examined by first comparing the mean drag profiles against reported wind-tunnel measurements that were carried out under similar test conditions and numerical simulations (LES and DNS). The surface drag is extrapolated from the nearest data point surrounding the roughness elements. Second, estimates of the friction velocity Uτp\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$U_{\tau }^p$$\end{document} and the zero-plane displacement height dp\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d^p$$\end{document} are obtained by integrating the axial pressure difference across each individual obstacle, assuming it is spanwise uniform. These are compared against direct measurements of the wall-shear stress from a floating-element balance and a pressure-tapped cube, as well as against estimates from indirect methods. In addition to mean pressure maps, snapshots of the pressure field are obtained via 2D-TH, based on Taylor’s Hypothesis, which are used to compute the RMS of the pressure fluctuations on the surface of a cube. The results indicate that 2D-RANS and 2D-TH perform adequately, providing reasonable estimates of the mean pressure distribution and of the boundary-layer flow parameters, outperforming indirect methods which rely on equilibrium assumptions that are often not verified.


Introduction
Experimental studies of boundary-layer flows over urban environments are generally limited to velocity data. Few provide a detailed description of the surface pressure over the roughness obstacles and it is not yet clear how direct measurements of the surrounding pressure field could possibly be achieved. In fact, only Cheng et al. (2007) and Claus et al. (2012b) have examined the change in pressure distribution and its effects on cube arrays with varying wind direction and packing density. Their analysis sought to document the flow statistics and the aerodynamic characteristics of different surfaces, i.e., the effective roughness height y 0 and the zero-plane displacement d, by pressure-tapping opposite faces of a roughness element. This approach is particularly attractive for d-type surfaces, like transverse ribs (Antonia and Luxton 1971), because a single wall-normal array of pressure taps is sufficient to estimate the form drag. In contrast, cubical arrays require a bidimensional grid, which, due to space restrictions, may have to be distributed over multiple identical obstacles to guarantee an adequate spatial resolution (e.g., Claus et al. 2012b). Further complications arise for surfaces with nonuniform height distributions, as each repeating unit is comprised of multiple elements with different aspect ratios, some of which may be too small to accommodate tapping holes. Measuring surface pressure then becomes a complex and somewhat cumbersome procedure. An alternative approach would be to use a non-intrusive technique to measure the surrounding static-pressure field and extrapolate the value on the surface. Direct methods have been successfully used, for example, using micro-air bubbles as pressure sensors (e.g., Ooi and Acosta 1984;Ran and Katz 1994), however, they are strictly used in water tunnels and still present serious limitations with respect to the density and the distribution of the bubbles. They are also not compatible with particle image velocimetry (PIV), so simultaneous measurements of the flow field and the overlying pressure distribution cannot be obtained. This would be of interest to further our understanding about the mechanisms responsible for drag generation and how these relate to spatial pressure fluctuations and large-scale motions above the canopy layer.
Given the lack of a suitable direct measurement method, pressure reconstruction from PIV stands as a promising tool. Introduced by Gurka et al. (1999), who estimated the average pressure fields of an impinging jet and of a channel flow, it quickly became popular for its potential to provide an integral description of the flow dynamics. At present, numerous variations of this technique have been proposed for planar, stereo or volumetric data, uncorrelated or time-resolved, using Eulerian (EU) or Lagrangian (LA) approaches. A review by Van Oudheusden (2013) details the operating principles and implementation aspects of these methods. He further discusses their intrinsic limitations arising from the seemingly chaotic motions of turbulent flows, and highlights the dominant factors that influence their accuracy, namely, the quality of the velocity-source data and their spatial and (if applicable) temporal resolutions. The Eulerian computation for convective flows using Taylor's hypothesis (TH), introduced by de Kat and Ganapathisubramani (2013), overcomes the requirement for time information to obtain estimates of the instantaneous pressure gradient. In this special case, performance is subject to the measurement of a convection velocity rather than the temporal resolution (Krogstad et al. 1998;Laskari et al. 2016).
Snapshot planar-PIV for uncorrelated in-plane velocity measurements has become widespread, but its potential for pressure reconstruction is essentially restricted to (quasi-)2D turbulence, when the contribution of out-ofplane motions to the transport of momentum is negligibly small. Yet, despite this detrimental aspect, some applications (e.g., industrial measurements) would largely benefit from full-field pressure estimations, even if the flow is not strictly two-dimensional. This premise has recently motivated the work of Van der Kindere et al. (2019), who quantified the accuracy of two PIV-based reconstruction methods using in-plane velocity data (2D-EU and 2D-TH), and how they perform against equivalent, time-resolved 3D approaches. They showed that, for a turbulent boundarylayer, using 2D-TH produces comparable results with those of a 3D-EU approach for higher noise levels in velocity data, and that pressure statistics are typically insensitive to the third spatial dimension. Furthermore, they achieved an agreement to within 4% ( ≈ 2 Pa) in the mean surface pressure over a forward-backward facing step with data from an array of surface-mounted transducers, and 10 − 15% in the root-mean-square (RMS) pressure. The correlation coefficient between the estimated and the reference pressure signals is approximately 0.5. In view of their analysis, employing these techniques towards understanding and modelling the flow processes in urban areas is a natural step.
In this study, static-pressure fields over two staggered arrays of cuboids are inferred from in-plane velocity data, using 2D-RANS and 2D-TH, by neglecting the contribution of the out-of-plane motions. This is especially detrimental within the canopy layer, where the flow is highly three-dimensional. Although empirical analyses suggest that the error induced by the missing terms is potentially the primary source of uncertainty, it cannot be readily quantified. The performance of the reconstruction methods is then examined following a type-B evaluation strategy: First, by comparing the mean surface pressure distribution against reported wind tunnel measurements (Claus et al. 2012b) and numerical simulations (Claus et al. 2012a;Leonardi and Castro 2010). Second, by estimating the values of friction velocity U and zero-plane displacement d based on the reconstructed pressure field. The pressure distribution is therefore assumed to be uniform along the spanwise direction, which, according to measurements of Cheng et al. (2007), is a reasonable approximation over nearly 70% of the width of the cuboids.

Experimental methods
Experiments were conducted in the open-return suction wind tunnel at the University of Southampton. This facility features a 7:1 contraction followed by a closed working section, 4 m long with cross section 0.9 m wide and 0.6 m high. Boundary layers are established directly on the wind tunnel floor, following the contraction. The freestream turbulence intensity has been reported to be homogeneous and lower than 0.5% at the measurement location (Claus et al. 2012b;Placidi and Ganapathisubramani 2018), and the flow develops in nominally zero-pressure gradient. The walls of the test section are optically transparent, to facilitate using imaging techniques (e.g., PIV), and the floor is cut out to fit a floating element (FE) balance. This work follows the convention that x, y, z are, respectively, the streamwise, wall-normal and spanwise directions; U, V, W, are the corresponding velocities in those directions.
We consider two rough surfaces, illustrated in Fig. 1: a staggered array of cubes with uniform height (C10U) and an array of cuboids with the same plan arrangement but a variable height distribution (C10R), the standard deviation is 3 mm (further details regarding its design can be found in Cheng and Castro 2002). Both arrays have the same mean height H = 10 mm and a planar solidity fraction p = 0.25 (defined as the ratio between the plan area of the obstacles and the floor area of a repeating unit). This particular kind of surface topology is often referred to as urban-like roughness because it shares similar features with urban environments. Particularly, the regular shape of the obstacles/buildings and the long street canyons, which, in this case, span across the wind tunnel. Cheng and Castro (2002), Reynolds and Castro (2008) and Claus et al. (2012b) have previously investigated the boundary-layer flow developing over these surfaces, under similar test conditions. They provide measurements of the vertical profiles of the mean streamwise velocity and of higher order statistics, obtained via cross-wire anemometry and Laser Doppler Velocimetry (LDV) . Reynolds and Castro (2008) additionally acquired snapshots of the velocity field over cube roughness using planar PIV. Since the present measurements were found to be entirely in agreement with these studies, only the reconstructed pressure fields are shown here, of which there is no previous available experimental data.

Floating element
Wall shear stress is directly measured using the FE balance described in Ferreira et al. (2018). Its design is based on the parallel-shift linkage, featuring dedicated pairs of bendingbeam transducers to monitor not only the streamwise load but also the induced pitching moment, as a means to decouple extraneous loads. The FE itself is a 200-mm-side square that is small enough in relative terms to ensure local measurements of wall shear stress. The balance is flush mounted with the wind tunnel floor, roughly 3.0 m downstream of the contraction, as depicted in Fig. 2. Special care was exercised mounting the surrounding tiles of the FE, given the tight tolerance of the fix-to-floating surface joint that is only   (3) is the sensing element. The field-of-view for the PIV is highlighted green and is located downstream of the balance = 0.5 mm wide. It is emphasised that the dimensions of the FE are not a multiple of a tile-size, which is 80-mm-side square. Therefore, on the floating element, two complete tiles were placed at the centreline, one behind the other, the first flush with the leading edge. Adjacent tiles were then laser cut to fit within its plan area while the remaining pieces were mounted outside, creating a seemingly continuous surface. Having an incomplete set of tiles on the sensing element does not impact the measurement of wall shear stress over the cube roughness, since there still is an integer number of repeating units. For the 'random' height case (C10R), however, the tile itself is the repeating unit so obstacles in rows 1, 2, 7 and 8, and columns C and D, are underrepresented (refer to Fig. 1). This potentially introduces some bias error that is challenging to estimate, but from the relative number of obstacles left out to the total number within 9 tiles covering the FE, we may assume their impact is not appreciable.
A total of five runs were performed in subsequent days for freestream velocities ranging between U 0 = 10 and 26 ms −1 in increments of 2 ms −1 . Each acquisition lasted 30 s at 150 Hz, equivalent to 2500 boundary-layer turnover times at the lowest operating speed. Pre-and post-calibrations were conducted for each configuration without notable discrepancies. The acquisition procedure and uncertainty analysis are detailed in Ferreira et al. (2018).

Particle image velocimetry
Measurements of the flow field that developed over each surface were acquired using planar PIV in the streamwise wall-normal plane. For the case of uniform height, the measurements were taken at 3 different spanwise locations (on the side of the cube, 1/4 cube and centre cube) to spatially average the flow across a repeating unit. Given the height variability of the second surface roughness, we considered instead 8 spanwise locations matching the centreline of the obstacles numbered according to Fig. 1. The Field of View (FOV) is located behind the FE, as shown in Fig. 2. This development length was verified to be sufficient for the flow to establish a constant shear-stress layer above the canopy and to reach fully rough conditions at a fixed freestream velocity of U 0 = 10 ms −1 ; Re | C10U = 5300 and Re | C10R = 6000 . Two laser/camera arrangements were setup depending on the surface roughness; Fig. 3 shows the schematic for the array of cuboids with a variable height distribution (C10R). The flow was seeded with ≈ 1 m particles of a vaporised solution of glycol-water. Two Litron 200 mJ dual pulse Nd:YAG lasers were used to illuminate the particles via a system of mirrors and 50 mm cylindrical lenses, which steer and expand the laser beams into sheets with a thickness no larger than 1 mm. The set of optics was mounted on a traverse that enabled accurate displacement of the light sheets in the spanwise direction to predetermined locations. The streamwise distance between them l was adjusted to eliminate the shadows cast by the obstacles in the canopy. Two LaVision Imager Pro LX 16MP cameras, each equipped with 200 mm f/8 Nikon lenses, allowed a FOV 0.8 wide and 1.2 high (in terms of mean roughness height, 9H × 14H ). Both cameras share the same FOV, but have different viewpoints to overcome the lack of optical access within the canopy. Due to space restrictions, one of them was fitted with a Scheimpflug adapter to adjust its plane of focus (the angle of the image plane relative to the plane of focus ≃ 10 • ).

Fig. 3
PIV setup. Lasers were placed above the test section. Two sets of optics and mirrors were mounted on a spanwise traverse, which allowed to simultaneously displace the light sheets to predetermined locations. The cameras were mounted on two-dimensional linear translation stages, outside the wind tunnel. One of them featured a Scheimpflug adapter to adjust the plane of focus to the image plane by ≃ 10 • . The streamwise distance between the laser beams l ensures the entire canopy region is illuminated Page 5 of 13 70 The setup of the cube roughness (C10U) employed a single camera and a laser in similar fashion as described above. For each test case, 2500 image pairs were recorded at a fixed frequency of 0.5 Hz, which is slow enough to safely assume statistical independence between them. The images were correlated in the pixel domain with LaVision Davis 8.2.2, using decreasing interrogation window sizes down to 16px×16 px and an overlap of 75% . The instantaneous vector fields were then filtered using a 3 × 3 Gaussian kernel and finally scaled and mapped to real dimensions, resulting in an effective spatial resolution of approximately 0.35 mm. The average free-stream pixel displacement is 16 and the particle diameter larger than 2p. Following Benedict and Gould (1996), the linear estimates of the uncertainty associated with the mean values are Velocity maps for C10R were obtained by combining two independent dual-pulse lasers whose light sheet overlap one another over most part of the FOV, except in the canopy. Variations in light intensity are not expected to introduce significant bias error, yet perfect alignment of the cameras and calibration of each individual FOV were paramount to stitch the vector fields (Raffel et al. 2018). The stitching algorithm uses a modified tapered-cosine window as the weighting function to ensure a smooth transition across the wall-normal boundaries of the overlapping region.

Pressure estimation
The numerical framework for pressure reconstruction methods has been outlined by numerous authors, so we only present here the general expressions for those employed in the current study, 2D-RANS and 2D-TH. The reader is referred to Van Oudheusden (2013) and Laskari et al. (2016) for an indepth explanation of the theoretical background, uncertainty analysis and implementation.
2D-RANS: Time-averaged pressure fields were estimated from PIV data using a Poisson formulation of the Reynoldsaveraged Navier-Stokes (RANS) equation, and assuming a divergence-free flow. Accordingly, where = + � is the velocity vector field, p is the average pressure field, R is the Reynolds stress tensor, is the density of the fluid and its kinematic viscosity. Equation 1 was spatially integrated by means of a Poisson solver, designed based on the work of de Kat and Van Oudheusden (2012) (velocity gradients are determined using a central difference scheme). Boundary conditions were enforced all around the domain of integration: along the top edge, which lies in the freestream, pressure values were prescribed by a modified Bernoulli equation to account for traces of turbulent kinetic energy. Around the roughness obstacles and on the vertical boundaries of the integration region, the pressure gradient was instead imposed via Neumann conditions. The freestream reference pressure p 0 was measured with a pitotstatic probe mounted at the same downstream location as the FOV.
2D-TH: A crucial aspect in determining instantaneous pressure fields from velocity information is the correct estimation of the material acceleration in the Navier-Stokes equation. For incompressible flows, and taking an Eulerian frame of reference, the latter may be written as follows Although most terms in Eq. 2 can be readily computed from snapshot PIV, the lack of temporal resolution impedes a complete description of the fluid motion. This is traditionally achieved with time-resolved measurements of the velocity field, using high-end imaging equipment or multiple-exposure techniques. Alternatively, de Kat and Ganapathisubramani (2013) proposed using Taylor's "frozen turbulence" hypothesis to approximate the motion of a perturbation (i.e., turbulent eddy) relative to the local mean as one of convection, thus relaxing the acquisition requirements. The time derivative of a velocity fluctuation then becomes where is the local convection velocity vector. After substitution into (2), this yields an expression independent of the time derivative.
The underlying assumption behind TH is that the time scale of a turbulent field is larger than the time scale of its advection downstream. This is typically true for convective flows, as in the case of grid-generated decaying turbulence, where turbulent eddies of different sizes are convected together with the local mean. In wall-bounded flows, however, turbulent structures have different characteristic transport velocities depending on the position within the boundary layer and on the scale and type of event (Krogstad et al. 1998), posing a practical challenge to accurately measure . Incidentally, researchers have often found the mean field to be a sensible estimate, yielding reasonable results at least for channel flows (see, e.g., Geng et al. 2015, Laskari et al. 2016) and forward-backward facing steps (Van der Kindere et al. 2019). Following this approach, = .

Implications of planar-PIV
The uncertainty associated with PIV-based pressure estimation is primarily related to the quality and nature of the velocity data. As explained by Van Oudheusden (2013), unless the flow is purely two-dimensional (2D), out-ofplane motion must be accounted for in the governing equations whether planar or volumetric pressure fields are to be extracted. In a statistical sense, the outer region of the boundary layer which develops over these surfaces may be treated as 2D, in which case neglecting 3D terms does not affect the accuracy of the estimation. Near the wall, however, the flow is dynamically influenced by the roughness length scales, giving rise to dispersive stresses. Especially within the canopy, where the flow meanders around the obstacles, strong out-of-plane motions are expected. There have been some attempts to circumvent the lack of 3D information in planar PIV (Haigermoser 2009;Baur 1999), but so far it appears that modifications to existing methodologies alone do not yield significant improvements (Van Oudheusden 2013). Therefore, additional effort would be necessary to resolve the missing components of velocity and acceleration (e.g., with Tomographic PIV or 3D particle tracking velocimetry).
Quantifying the error introduced by this approximation is nontrivial, as it requires a priori knowledge of the flow conditions. Sensitivity studies (Charonko et al. 2010;de Kat and Van Oudheusden 2012) have shown, nonetheless, that a small to moderate degree of out-of-plane motion does not seriously affect the pressure field estimation, provided a global integration approach such as a Poisson solver is used-the error within the canopy is diffused into the outer region where the flow is effectively 2D. The impact of 3D flow motion on planar reconstruction was further investigated by McClure and Yarusevych (2017). They found that for fully turbulent regimes, the relative error in the near wake of a cylinder does not exceed 5% , but can grow larger than 20% farther downstream, as turbulence becomes increasingly homogeneous (and in-plane vorticity emerges). On the basis of these and previously existing empirical analysis (e.g., de Kat and Van Oudheusden 2012), they concluded that planar methods are generally applicable if the out-ofplane gradients are less than half of in-plane velocity gradients, from which point the accuracy of the estimation rapidly decays. To assess this criteria, we use the large eddy simulation (LES) data set of Xie and Castro (2006) to compute the relative magnitude of the mean velocity gradients at the centreplane of a cube, shown in Fig. 4 for the streamwise component. It is clear that U∕ z is consistently smaller than U∕ x or U∕ y by one order of magnitude, occasionally exceeding the aforementioned threshold. The uncertainty associated is then expected to remain at acceptable levels in the region of interest. Naturally, above the canopy the streamwise gradient approaches zero causing the ratio to suddenly diverge.

Additional error sources and noise propagation
Besides the lack of spatial information, 2D-TH formulation is intrinsically sensitive in regions of intense shear and turbulence intensity. This is partly due to the limited spatial resolution of the PIV system but, more importantly, because the basic assumptions of TH are unlikely to hold. This limitation was recently investigated by Van der Kindere et al. (2019) for a turbulent-boundary layer flow past a wall-mounted rib. TH appeared to perform reasonably well over most part of the FOV, except in the vicinity of the shear layer shed off from the leading edge and in the wake region, where the mean velocity is considerably lower than the local fluctuations. In light of (a) (b) Fig. 4 The missing non-linear convective terms. The mean ratio between in-plane and out-of-plane streamwise velocity gradients within the canopy of a staggered-cube array, from LES data of Xie and Castro (2006). Slice along the centreplane of the cubes. a = U∕ x and b = U∕ y . Gray colormap is for values larger than 0.5 these result, they suggested | |∕| | ⋅ H∕U 0 > 10 to be a suitable performance indicator for 2D-TH, where is the mean out-of-plane vorticity and H the obstacle height. Accordingly, using the mean roughness height as the length scale to normalise vorticity, we expect TH to yield sensible results above the canopy. Values are typically less than two in the shear layer and below the canopy top, except in the vicinity of isolines of zero-velocity where it far exceeds the stipulated limit.
Noise propagation from velocity data was estimated following the procedure outlined by de Kat and Van Oudheusden (2012), which is derived from the approach to uncertainty analysis of Kline and McClintock (1953). For 2D-RANS, neglecting the contribution of the viscous term, the linear propagation of uncertainty applied to Eq. 1 yields where is the estimated precision error of the arbitrary quantity and xy is the spatial resolution of the vector fields, listed in Table 1. Note that Eq. 5 considers secondorder central finite differences to estimate derivatives in space; truncation errors were not included in this analysis. For 2D-TH, Laskari et al. (2016) expressed the relative error propagation from velocity as We estimate 2 p RANS to be on the order of 2% around the shear-layer interface and in the region immediately above, where turbulence intensity is highest (refer to Cheng and Castro 2002). Within the canopy, it becomes noticeably smaller (less than 1% ). Lastly, to determine p TH , we assume the uncertainty on the local convection velocity u c ≈ 3% , yielding an average relative error in pressure of about 1.2%.

Analysis of the pressure fields
Mean pressure maps, reconstructed from velocity data using 2D-RANS are shown in Figs. 5a and 6, expressed nondimensionally by the pressure coefficient C p = (p − p 0 )∕q 0 , where q 0 is the freestream dynamic pressure; the flow direction is from left to right. Within the canopy, the pressure field is governed by the presence of roughness obstacles. High pressure regions generally develop on the upper half of the windward side, where the mean streamwise velocity is highest. On top, a low pressure region reveals the existence of a shear layer shed off the sharp leading edge that is essentially stronger for tall isolated obstacles. The pressure (a) (b) Fig. 5 Pressure within the canopy of C10U. a The normalised mean pressure field reconstructed from planar-PIV on the vertical centreplane of the obstacles using 2D-RANS, and b the absolute difference between the estimates obtained with 2D-RANS and 2D-TH.
The pressure coefficient C p = (p − p 0 )∕q 0 , where q 0 is the freestream dynamic pressure. Obstacles in the measurement plane are highlighted gray is quickly recovered downstream, suggesting that the total force acting on the cuboids is predominantly influenced by the pressure distribution on the windward side. Effects of mutual sheltering are particularly visible over C10R (Fig. 6), whereby the pressure surrounding the smallest obstacles is markedly modified. Above the canopy, the pressure field is mostly uniform and takes on the value of the local freestream static pressure p 0 ≃ 0. The absolute difference between estimates obtained with 2D-RANS and 2D-TH, averaging over 2500 snapshots of the pressure field, is given in Fig. 5b for C10U. Discrepancies are most significant around the bottom-left corner, reaching 20% of the largest pressure value, but appear to be confined to a small region of the FOV. They presumably emerge from noise propagation and the sensitivity of the Poisson solver to Neumann boundary conditions. Along the edge surrounding the obstacles, which is typically 0.5 mm off the wall, both methods produce similar results. In fact, the spatially averaged absolute difference ⟨�C p,RANS − C p,TH �⟩ is only 0.0038, comparable to values reported by Van der Kindere et al. (2019). Despite having been inferred from distinct PIV data sets, the pressure fields over C10R agree at least qualitatively with each other, in the sense that the higher pressure regions pertain to the tallest obstacles (Fig. 6). They are also not significantly affected by the mask of spurious vectors (in white), except for (d) and (g), where there seems to exist . 6 Pressure within the canopy of C10R. Normalised mean pressure fields obtained using 2D-RANS at each streamwise alignment, 1-8, as outlined in Fig. 1. x-coordinate is referenced to the leading edge of the roughness tile. Obstacles in the measurement plane are highlighted gray a local artificial decrease of the pressure value in-between the mask of the out-of-plane obstacles (at x∕H = 3.5 and 4.5, respectively).

Drag profiles and pressure fluctuation
To assess the performance of the pressure estimation, we examine the normalised axial-pressure difference across selected obstacles p(y) = p w (y) − p l (y) , where the subscript identifies the windward (w) and leeward (l) faces, represented in Figs. 7a and 8. Available experimental measurements (Cheng and Castro 2002) and numerical solutions (Claus et al. 2012a;Xie et al. 2008;Leonardi and Castro 2010) were included for reference, revealing a notable collapse in spite of the inherent uncertainty of the reconstruction methods. Around the mid-canopy height of C10U, the agreement is generally better than 0.1⟨ p⟩ , but it worsens near the wall. There, a low pressure level is maintained by the recirculating region that develops ahead of the cube, so p is relatively small. It then grows larger away from the wall, reaching a peak for y > 0.8h . Leonardi and Castro (2010) showed that, for staggered cube arrays, the peakto-peak amplitude increases with packing density, as the shear-layer interface between the canopy and the flow aloft becomes gradually more prominent. Presumably, a similar effect would result from an increase of the roughness Reynolds number based on the mean height H + = U H∕ , which is somewhat lower for computations than it is in this experiment. Since the peak-to-peak amplitude is smaller with 2D-RANS, discrepancies in local extrema likely arise from uncertainties, both in the measurements and numerical solutions. C p rms profiles are shown in Fig. 7b. These are naturally prone to a higher noise propagation and are presented here without quantitative validation, though results from Van der Kindere et al. (2019) indicate that the relative standard error should lie within 10 − 15% . Pressure fluctuations are most intense around the stagnation region, on the windward side close to the canopy top, reaching nearly 50%  Fig. 1) are, respectively, 1.71H, 1.36H and 1H high. Colors identify different streamwise alignments Form drag 0.560 -Clauser 0.490 − 0.79 -of the local mean pressure value, and become less significant towards the wall. On the leeward side, C p rms is mostly uniform and relatively small, which reflects the existence of a low momentum region that extends across the canopy height.
Estimates of drag profiles for three cuboids within C10R are compared in Fig. 8a against LES data from Xie et al. (2008). We emphasise that the simulation was performed over a geometrically identical array, but the flow is in the opposite direction. Despite this important caveat, their results are still consistent with the current measurements, suggesting that the pressure field is hardly affected by the local height distribution. This is true, at least, for staggered and aligned arrays, provided the plan solidity fraction p ⩽ 0.25 , conditions which minimise transverse wake interaction. The agreement between drag profiles is best for the shortest obstacle (8-A). So it is reasonable to assume that effects of local Reynolds number become noticeable for taller obstacles, explaining, albeit not entirely, the reduced peak-to-peak amplitude from LES data (carried out at a lower H + ).

Boundary-layer flow parameters
Neglecting the viscous contribution to wall-shear stress, we attempt to estimate the friction velocity U over C10U and C10R from pressure data. The results are compared against FE measurements, listed in Table 2, to ascertain the performance of the reconstruction methods in the canopy layer. Additional estimates were included for reference, specifically those inferred from cross-wire anemometry of the boundary-layer profile by Cheng and Castro (2002), over identical obstacle arrays, and from a pressure-tapped cube (C10U) by Claus et al. (2012b).
The form drag of individual obstacles F i is obtained by integrating the axial-pressure difference over their cross sectional area h i × w , assuming a spanwise-uniform distribution and a zeroth-order (constant) extrapolation of the pressure value at the edge of the domain down to the wall. The surface shear stress is then estimated considering the total drag within the plan area of a repeating unit S. Accordingly, and The results, listed in Table 3, are consistent with direct measurements of wall-shear stress using the FE balance ( Table 2). The relative discrepancy between them does not exceed 3.7% with 2D-RANS and 7% with 2D-TH (C10U). Marginally lower pressure-based values are expected as the frictional drag is not accounted for, yet the inherent uncertainty level frustrates any attempt to quantify its contribution. Especially considering that U p relies on information at the vertical centreplane of the cuboids alone and, thereby, is likely overestimated. The total stress method appears to systematically under-predict the friction velocity, which could either result from the lack of a fully-developed equilibrium layer or, as Womack et al. (2019) recently argued, due to a small favourable pressure gradient imposed by the fixed cross-section of the wind tunnel facility. The zero-plane displacement d could also be determined following the theoretical arguments of Jackson (1981), who proposed that the vertical origin of the boundary layer is given by the centroid of the distributed drag. This reads, For the uniform array, values of zero-plane displacement agree with pressure tap measurements by Cheng et al. (2007) ( d p = 0.612H ) and Claus et al. (2012b) ( d p = 0.56H ), as well as with DNS data by Leonardi and Castro (2010) using pressure-based ( d p = 0.617H ) and log-law fitting methods ( d = 0.645H ). Note that Cheng et al. (2007) also present estimates obtained via a modified Clauser chart using U inferred from the total-stress method, all of which are greater than 0.822H. Additional estimates by Claus et al. (2012b) lie in the range 0.49H − 0.79H , depending on the slope of the log-law-they compared multiple values of U and, in some cases, treated the von Kármán coefficient as a fitting parameter. The marked discrepancy between reported values indicates that the uncertainty associated with indirect methods is in fact larger than is usually reported, impeding for example a definite conclusion on the validity of Jackson's Hypothesis.
The influence of height variability on surface roughness was systematically investigated by Jiang et al. (2008), Hagishima et al. (2009), andMillward-Hopkins et al. (2013), invariably leading to an increase of wall-shear stress and (7)

Distribution of surface drag over C10R
Previous studies (Kanda 2006;Xie et al. 2008;Millward-Hopkins et al. 2013) have shown that the distribution of the surface drag is dominated by the height variability of the roughness elements. Particularly over C10R, Xie et al. (2008) predicted that the percentage contribution of the tallest obstacle alone is 22.4% , despite its proportionate cross sectional area ( 10.8% ), whereas the second tallest are responsible for 42.9% and the remaining obstacles for 34.7% . An equivalent analysis can be done using Eq. 7 to estimate the relative contribution to the total drag of individual roughness elements. Illustrated in Fig. 9a, the results show a remarkable quantitative agreement with findings of Xie et al. (2008). The load exerted by the fluid is fairly uniform for obstacles which stand above the mean canopy height. 2-A, 7-B, 6-C and 5-D protrude 0.36H and produce about the same form drag, in the range 9.5 − 11.3% . In contrast, the variability is increased for those which lie in sheltered regions. Depending on the local height distribution, 1H-high obstacles experience loads from 0.9 (3-D) up to 9.5% (1-B) of the total surface drag.

Concluding remarks
We investigated the potential of using two-dimensional, PIVbased pressure reconstruction methods to achieve a more complete description of the flow field over two staggered arrays of cuboids (C10U and C10R). Empirical analyses suggest that the error induced by the missing out-of-plane components of velocity and acceleration could be problematic within the canopy layer, where the flow meanders around the roughness obstacles, however, it cannot be quantified without a priori knowledge of the flow conditions. So to assess the performance of this approach, surface pressure was extrapolated from the edge of the domain and compared against equivalent numerical and experimental studies. Overall, both 2D-RANS and 2D-TH yielded sensible estimates over the entire domain, showing minor discrepancies that should arise due to the sensitivity of the Poisson solver to boundary conditions. Estimates of the mean axialpressure difference across individual elements showed a remarkable collapse with LES and DNS data, as well as with pressure-tap measurements. The RMS of the pressure fluctuations on the surface of a cube (C10U) was also obtained via 2D-TH, although it could not be validated against published results. This may potentially be explored to study, in a statistical sense, the relationship between the streamwise load applied on the roughness obstacles and the turbulent structure of the boundary layer above. It would then be beneficial to quantify the correlation coefficient between pressure values from PIV and reference transducer data on the a Plan-view, as given in Fig. 1, wind direction is from left to right. Each individual obstacle is color-coded according to the relative contribution to the total surface drag, obtained by integrating the axial pressure difference over the cross area. Individual percentage contributions are written either above or below. b Pie chart illustrates the relative contribution to wall drag by the tallest obstacle (1.72H high) and by the second tallest (1.36H high) surface. Lastly, the frictional velocity U and the zero-plane displacement height d were inferred assuming a constant spanwise distribution of the surface pressure. Estimates closely match direct measurements using a FE balance and a pressure-tapped roughness obstacle. Discrepancies are at worst 7% in U and 3.5% in d.
We conclude that pressure reconstruction from planar velocity data, using 2D-RANS and 2D-TH, has the potential to be a very useful tool in the study of urban boundary layers. It revealed a complex interaction between roughness elements, having different drag profiles, depending mostly on the relative height distribution along the streamwise direction. This is reflected on the contribution to wall drag by individual elements that is largely disproportionate over C10R, some of which might even be thrust producing. This approach additionally stands as an alternative to direct measurement techniques when these are not available or cannot be easily employed (such is the case of C10R, whose repeating unit is comprised of 16 different elements), outperforming an ill-conditioned three parameter fit of the mean velocity profile or the total stress method. There are, nonetheless, some rather limiting aspects to bear in mind: For arbitrary-shaped obstacles (e.g., pyramids, cylinders, spheres) or arrays of cuboids at varying wind directions, the contribution of the spanwise components is likely to become important, resulting in a less accurate estimation in the canopy layer. It also requires highly spatially resolved velocity maps and, since surface pressure is extrapolated from the closest available data point, performance is subject to laser reflections off the wall. Some of these shortcomings may be mitigated, however, using thin-volume tomographic PIV or a refractive-index matched facility.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.