Universal scaling-law for flow resistance over canopies with complex morphology

Flow resistance caused by vegetation is a key parameter to properly assess flood management and river restoration. However, quantifying the friction factor or any of its alternative metrics, e.g. the drag coefficient, in canopies with complex geometry has proven elusive. We explore the effect of canopy morphology on vegetated channels flow structure and resistance by treating the canopy as a porous medium characterized by an effective permeability, a property that describes the ease with which water can flow through the canopy layer. We employ a two-domain model for flow over and within the canopy, which couples the log-law in the free layer to the Darcy-Brinkman equation in the vegetated layer. We validate the model analytical solutions for the average velocity profile within and above the canopy, the volumetric discharge and the friction factor against data collected across a wide range of canopy morphologies encountered in riverine systems. Results indicate agreement between model predictions and data for both simple and complex plant morphologies. For low submergence canopies, we find a universal scaling law that relates friction factor with canopy permeability and a rescaled bulk Reynolds number. This provides a valuable tool to assess habitats sustainability associated with hydro-dynamical conditions.

In the last decades, an increasing number of countries have introduced restoration practices to promote the sustainable management of rivers 1 . These practices often advocate the use of in-channel and riparian vegetation to preserve essential ecosystem services such as the reduction of banks erosion and the sheltering of fishes and invertebrates 2 . One key factor to properly assess water quality and ecological functions of rivers is the quantification of flow resistance induced by vegetation [3][4][5] . By creating an additional drag, canopies reduce the flow velocity in the vegetated layer, promote local retention of chemicals and deposition of particles that improve the stability of riverbanks and affect the fate and transport of nutrients and pollutants along the river [6][7][8] . While momentum and mass transfer mechanisms between the canopy and the free flows chiefly depend on vegetation geometry [9][10][11][12][13][14] , the characterization of friction in terms of vegetation morphology and flow regimes is an open challenge. Many studies have investigated the relationship of the drag coefficient C D with Reynolds number 5,[15][16][17][18][19][20][21] , stem population density 15,22,23 , canopy spacing 15,21 , plants geometrical arrangement 24,25 , blade flexibility 25,26 , leaf area to stem area ratio 26 , solid volume fraction 21,27,28 , penetration length [29][30][31] and the Keuglan-Carpenter number for wave-dominated flows 16,32,33 . Yet, universal parametrizations of the drag force through a drag coefficient C D have proven elusive due to its dependence on many parameters, notably the combination between canopy morphology and dynamic conditions. In alternative, empirical coefficients, e.g. Manning friction factor, can be used to estimate the flow resistance. However, their applicability is generally limited to the conditions in which they are originally estimated and their predictive capability has often been criticized, especially when dealing with complex channel and plant morphologies 34 . Natural vegetation can exhibit both relatively simple geometries (e.g. in grasses) and complex shapes as in plants with foliage, branches and leaves that create a vertical variability in the plant density 12 . This renders quantification of the compound effects of different plant components on flow resistance even more dire 20,[35][36][37] . Leaves and branches increase the total surface area and the drag force 2,4,5 , while this effect decreases for flexible vegetation due to the reconfiguration of the plant in the streamline direction 38,39 . As a result, the relationship between the drag force and the mean flow velocity is in general linear or less than quadratic (with an exponent between 1.3 and 1.9) 4,19,40 .
A variety of flume studies with real plants or plant prototypes mimicking natural vegetation have been conducted to investigate the effects of natural vegetation on flow 38,39 . These studies suggest that the most important

Model Description
Following the modeling framework described by Battiato and Rubol 11 , we consider a two-dimensional fully developed incompressible turbulent flow in an open channel of slope θ with S 0 := tanθ ≈ sinθ, whose bottom part, ∈ z H (0, ), is occupied by a morphologically complex canopy layer of (deflected) height H. The flow depth is H + L. A sketch of the model is shown in Fig. 1. The canopy is treated as a porous medium of permeability K [L 2 ]. The effective permeability K, an intrinsic property of an obstruction and its morphology (i.e., porosity, density etc.), is a metric that describes how easily water flows through the canopy. The higher the permeability the lower is the resistance of the porous medium to flow. Similarly, the lower the permeability, the greater is the resistance to flow offered by the vegetation. In the context of classical porous media (e.g. consolidated and unconsolidated granular media like rocks and soils), there are numerous empirical relationships that allow one to estimate K from, e.g., porosity. One such famous relationship is the Karman-Cozeny equation. Yet, one should be cautious about the use of such empirical formulas to canopy layers since the structure of the obstruction is significantly different from that of geologic porous media (from spheroidal obstacles to fibrous-like obstacles) 51 . We emphasize Figure 1. Sketch of the conceptual model. The deflected vegetated layer, of height H, is characterized by an effective permeability K, which describes the ease with which water can flow through the canopy layer. The water level above the canopy height is L, and H + L is the flow depth. The mean velocity profile above and within the canopy is ˆû z ( ). The triangle indicates the location of the free surface, i.e. the air-water surface where the shear stress is zero. that the link between canopy permeability K and geometrical features of the vegetated layer (e.g. canopy density, porosity, leaf area density etc.) has been explored by Battiato and Rubol 11 for canopies constituted by regular arrays of rigid cylinders. The extension to flexible canopies with complex morphology is object of current investigations. For completeness, we report the model equations developed by Battiato and Rubol 11 . The mean velocity profile above and within the canopy layer is obtained by coupling the log-laŵˆˆˆˆκ to the Darcy-Brinkman equation , and κ is the reduced von Kármán constant 52,53 (see 11 for details). We define the following dimensionless quantities where q = ρgS 0 H 2 /μ e is a characteristic velocity scale. The analytical solution of Eqs (1-2) for the mean velocity profile û z ( ) inside and above the canopy layer is 11 respectively, where The Darcy friction factor f is defined as Among the relevant parameters that are used to quantify the flow resistance, we report the expression for the inverse of the drag length scale C D a, given by the product of the medium drag coefficient, C D , and the frontal area of the canopies, a Once the channel morphology (i.e. H, L, B, S 0 ) and the canopy effective properties (i.e. its permeability K) are known, Eqs (4), (6) and (7) can be used to analytically calculate/predict the mean velocity profile above and within the canopy û z ( ), the volumetric discharge Q and the Darcy friction factor f. In the following Section, we compare model predictions and experimental data published in the literature.

Results and Discussion
In this Section, first we present a comparison between model predictions and experimental measurements collected from the literature, and then identify a universal scaling law for the friction factor across different canopy morphologies. The data include velocity profiles, flow rates and friction factors collected both in flumes and real channels, and span a variety of vegetation types and morphologies ranging from grasses (e.g., eelgrass, cordgrass and barley) to woody vegetation (e.g., poplars) and bushes (e.g., willows), with plant heights varying from 0.05 m Velocity profiles. To test the capability of the proposed model to reproduce the flow within and above canopies with different morphologies and flexibilities, we first compare the analytical solution in Eq. (4) to the experimental datasets listed in Table 1. Figures 2 and 3 show the comparison between the measured and predicted velocity profiles. The comparison is conducted by fitting either one (i.e. canopy permeability, K) or two (i.e. canopy permeability, K, and the reduced von Kármán constant, κ) parameters using the orthogonal distance algorithm (ODA 54 ), which minimizes the orthogonal distance between the model curve and the experimental data (solid black and grey lines for the one and two parameter fitting, respectively). Figures 2 and 3 show that the model (both with 1 and 2 fitting parameters) successfully reproduces the mean velocity profiles in the canopy layer across a wide range of canopy morphologies for both real plants and plastic prototypes. Additionally, surface layer velocity predictions with the one-fitting parameter model improve with increasing submergence, as shown in the profiles of Nepf and Vivoni 55 (NV125, NV150, NV175, NV190 and NV275) in Fig. 2 and Cassan, et al. 56 (CassanS11, CassanS12, CassanS21, CassanS22) in Fig. 3. Remarkably, despite the model is developed under the hypothesis of non-deformable media, it is still capable to well reproduce the velocity profiles measured in flexible vegetation, including tall willows, i.e. RSparse and RDense in Fig. 2 (dataset of Righetti, et al. 57 ), in patches of poplars, i.e. SiniscalchiH in Fig. 3 (dataset of Siniscalchi, et al. 58 ), as well as the velocity profiles collected by Cassan, et al. 56 in a real channel with prone non-uniform vegetation, i.e. CassanS11, CassanS12, CassanS21, CassanS22 in Fig. 3.
While the one-fitting parameter model always successfully captures the flow field in the canopy layer for the whole dataset of Figs 2 and 3, the two-fitting parameter model is necessary to improve the description of the velocity field in free-surface layer in channels where the water height above the canopy constrains the development of the log-law (see the profiles RDense of Righetti 57 , NV125 of Nepf, et al. 55 and VELASCO of Velasco 59 in Fig. 2). This is expected since for low submergence the boundary effects at the surface layer become important throughout the flow domain, and an additional fitting parameter (κ) is needed to compensate for the lack of accuracy of the log law at the free surface. The average value of the von Kármán constant in the two-fitting parameter model is κ = 0.22 (for the low submergence cases), only slightly higher than the value κ = 0.19 used for the one-fitting parameter model. While an open debate still exists on the universality of the von Kármán constant, we find that the fitted values are consistent with the values measured in vegetated flows 53 (see 11 for additional discussion). Instead, one fitting parameter is sufficient when the water height above the canopy is larger than the canopy  56 (i.e. CassanS11, CassanS12, CassanS21, and CassanS22) refer to highly submerged canopies (i.e., submergence greater than ten) for which the velocity profile is characterized by boundary layer rather than mixing layer features 2 . The model well reproduces the shape of the measured profiles with only slight deviations for points in proximity of the free surface, for which  Table 1. In the Figure, u is the dimensionless mean velocity, z is the dimensionless coordinate. Their definition is provided in Eq. (3). The solid black and grey lines correspond to the dimensionless mean velocity profile u predicted by (4) using one and and two fitting parameters, respectively. The list of parameters used for the fitting is provided in Table 2. The insets contain sketches of the canopy used in each experimental study: color and black drawings indicate real and plastic vegetation, respectively. presence of secondary currents causes a maximum of velocity below the air-water interface (see Cassan et al. 56 for details). Similarly, deviations between the data and the two-fitting parameter model were observed close to the free surface for the Shi40 dataset 60 in Fig. 3 and the NV275 dataset 55 in Fig. 2, where the experimental profiles capture the no-stress condition on the free surface 2 . We emphasize that these deviations between data and model are expected since the log law is not accurate at the air-water interface. While the model is able to capture the velocity profiles of flexible vegetation considered in this study, it cannot, in its current form, predict secondary  Table 1. In the Figure, u is the dimensionless mean velocity, z is the dimensionless coordinate. Their definition is provided in Eq. (3). The solid black and grey lines correspond to the dimensionless mean velocity profile u predicted by (4) using one and and two fitting parameters, respectively. The list of parameters used for the fitting is provided in Table 2. The insets contain sketches of the canopy used in each experimental study: color and black drawings indicate real and plastic vegetation, respectively.
ScIentIfIc REPORTS | (2018) 8:4430 | DOI:10.1038/s41598-018-22346-1 maxima of velocity due to strong vertical variability in plant density induced, e.g., by reconfiguration. Nevertheless, the model can be generalized to account for vertical heterogeneity and poro-elasticity by (i) assuming that permeability varies along the vertical direction (e.g., K = K(z)) and (ii) coupling the plant deflection with the flow field, following an approach similar to that proposed by Battiato, et al. 47 , who modelled the bending of carbon nanotubes forests under aerodynamic shear.
The fitted values of permeability K vary in the range 10 −4 -10 −1 m 2 or between 0.6-5 in term of the inverse dimensionless permeability, λ = H K : / . As expected, low permeable porous media are associated to highly dense canopy as for the dataset VELASCO 59 in Fig. 2, while highly permeable media correspond to the sparse willows of Righetti 57 (RSparse and RDense in Fig. 2). Shucksmith, et al. 61 (datasets Shuck7 and Shuck10 of Fig. 3) measured the velocity profiles at different stages of the plant growth and found that the velocity inside the channel decreased as the plant height and the plant density increased. This is also consistent with the increase in friction values described in the following section. A similar effect was also found for the data of Shi 60 (Shi20, Shi40 and Shi60 in Fig. 3). Changes in permeability due to changes in plant density were observed also for the dataset RSparse and RDense of Righetti 57 , for which the permeability values more than halved from the dense (2.5 bushes per square meter for the RDense dataset) to the sparse (1.1 bushes per square meter for the Rsparse dataset) configuration of willows. A summary of the relevant data including the value of the fitting parameters and the error between the analytical model and the experiments is provided in Table 2.
To further test the predictive capabilities of Eq. (4), we use the model calibrated on six datasets from Figs 2 and 3 (NV150, WF, VELASCO, BR2, SiniscalchiH and CassanS21) to perform a series of fit-free predictions on a different set of data (NV175, WF2, VELASCO2, BR3 and SiniscalchiM, CassanS22, respectively), i.e. we use calibrated permeability and von Kármán constant values to perform fit-free predictions of the mean velocity profiles collected for plants with the same morphology but under different dynamical conditions (i.e. water level and volumetric discharge). Figure 4 illustrates the comparison between the measured and the predicted mean velocity profiles, and shows good agreement within the uncertainty level (10% of the permeability value) for both morphologically simple (e.g., NV175, BR3) and complex (e.g., WF2, SiniscalchiM) canopies. Predictions of the velocity profiles within canopies with simple morphologies (e.g., BR3 62 ) were more accurate than those for   Table 2 provides the parameters values resulting from fitting the canopy permeability K only. The right portion of the Table 2  morphologically complex canopies (e.g., SiniscalchiM 58 and WF2 63 ). This is expected because the projected area (and therefore, the permeability) of plants with complex morphology is more sensitive to changes in water level than plants with simpler shape. This is due to the presence of leaves and branches that become more streamlined as velocity increases 64-66 . Discharge and Friction Factor. Once the fitted (Figs 2 and 3) and predicted (Fig. 4) velocity profiles are determined, we use Eq. (6) to calculate the discharge. Figure 5(left) shows the comparison between the measured and predicted discharge, Q exp and Q cal , respectively. The grey symbols refer to the purely predicted discharge (i.e. no fitting parameters) obtained from the profiles of Fig. 4, while the remaining symbols have been obtained with the one-fitting parameter model (i.e. solid black lines in Figs 2 and 3). The good agreement between data and predictions demonstrates that the one-fitting parameter model is sufficient to predict the flow rate in vegetated channels with sufficient accuracy. The two-parameter model can be used when highly precise measurements of the flow profile are required.
The vegetative resistance, i.e. the Darcy friction factor f (=C f /4), can be estimated from Eq. (7). Figure 5(right) shows a good agreement between the experimental and analytically predicted Darcy friction factors, f exp and f cal , respectively. We emphasize that, in this formulation, the information about the canopy morphological complexity is entirely accounted for through the effective canopy permeability K (or its dimensionless counterpart λ −1 ). As a result, f depends on canopy morphology through K and its direct effect on modulating the flow rate Q in the channel as evident from Eq. (6). Maximal values of friction factor are obtained for the plastic prototypes and for the highly dense (22500 stem/m 2 ) canopy of Velasco et al. 59 (VELASCO and VELASCO2). Among the prototypes, the  Cassan 56 (CassanS11, CassanS12, CassanS21, CassanS22). This is in line with the experimental work of Velasco 67 , that reported minimum resistance to the flow for highly deflected vegetation, with a friction factor comparable to the value of non-vegetated channel, f = 0.15, measured in gravel beds. The above considerations are also consistent with the medium drag C D a values calculated from (9). Noticeably, the computed values fall within the same order of magnitude of the values reported by Vargas 42 (see their Table 4), with higher medium drag for grasses than shrubs and forests. Also, the increase of shear stress with increasing plant height observed in the dataset of Shi 60 (Shi20, Shi40 and Shi60) is in agreement with the increase of the Darcy friction factors.

Universal Scaling of the Friction Factor. A number of works have pointed out that dynamic similarity
can be often observed in systems characterized by coupled flows through and over permeable media 11,48,68,69 . After model validation, we use the porous medium analogy to investigate the existence of such a universal scaling law for vegetated flows across different canopies morphologies. Our main result is the derivation, detailed in the Methods section, of the universal scaling between the Darcy friction factor, Reynolds number and canopy permeability.
The drag coefficient is defined as where D is the drag force and  is a characteristic length scale of the system. Using Bunckingham's Pi theorem and postulating the existence of self-similarity of the second kind (see Methods Section), we find that C D scales a D when the canopy layer is thick (or for shallow submergence), i.e. δλ Λ = = ≤ L K : / 1 (see Methods Section for detailed derivation). In (11), , we plot f/λ as a function of Re  where f is the experimentally measured friction factor, α = 1 and  = + H L, see Fig. 6(left). While it is expected that the global resistance in vegetated channels decreases with the  Re number, as plants become more streamlined to the flow, Fig. 6(left) shows that the friction factor exhibits the universal scaling proposed in (11) across different canopy morphology. Specifically, f/λ scales with the inverse of Re  , the modified Reynolds number, for thick canopies (i.e. Λ ≤ 1, black points in Fig. 6(left)).
The data regression of f/λ measured under shallow submergence (Λ ≤ 1) results in the following expression The proposed scaling (13) suggests that the friction factor appropriately normalized by the dimensionless permeability follows a universal scaling for thick canopies (Λ ≤ 1, or canopies with small submergence). The deviation of Cassan data (grey points in Fig. 6(left)) is expected since they fall in the category of thin porous media (Λ  1) or canopy with high submergence, where boundary layer-type velocity profiles develop. The fitted exponent β ≈ 0.86 corresponds to a less than quadratic relationship between the drag force and the mean flow, specifically ∼ .
. This is consistent with the scaling experimentally measured by Sand-Jensen 40 , who found exponents in the range (1.3, 1.9) and Armanini et al. 19 , who found an almost linear relationship between the drag force and the velocity for flexible submerged canopies.
Since  λ = (1) for the data analyzed, we hypothesize that λ = (1)  for most submerged canopies. Under this hypothesis, (13) can be further simplified to In Fig. 6(right) we plot the experimentally measured friction factor f in terms of Re  for a much wider dataset including rigid and flexible vegetation. Specifically, Fig. 6(right) includes the dataset of Table 1 in addition to 132 data points analyzed by Poggi, et al. in their Table 1 Table 1: f is the Darcy friction factor, defined in Eq. (7), λ = H K : / is the dimensionless inverse canopy permeability and Re  is the rescaled Reynolds number defined in Eq. (14). The black and grey symbols correspond to thick (low submergence) and thin (high submergence) canopies, respectively. The inset shows the same comparison is a semilog scale.

Conclusions
Accurate estimates of the impact of vegetation on flow resistance in riverine systems are critical to assess river management and restoration and to ensure aquatic habitat resilience and health. A common way to evaluate the quality of aquatic habitats is to determine whether selected hydraulic parameters (such as the water depth and flow velocity) fall within a range considered suitable for the biological needs of a specific organism. Since submerged vegetation reduces the mean velocity in the channel, accurate predictions of vegetation-mediated flow resistance is crucial to monitor the biotic diversity, to preserve a healthy life for the macro-invertebrates and fishes that find shelter in the submerged canopies, and to control the deposition of inorganic sediments in marshes 71 . Characterization of flow resistance for channels with topologically complex vegetation has proven elusive, though. On one hand, attempts to correlate the drag coefficient to relevant flow and vegetation parameters has opened a Pandora's box in which the intertwining of topological and dynamical controls on friction has been hard to disentangle. Similarly, simple experimental formulas such as Manning's, have proven to have poor predictive capability, especially in channels with complex morphology 34 . On the other hand, a few studies have demonstrated unique universal features of flows over obstructions of disparate types and topologies 48,68 ranging from submerged aquatic vegetation canopies, terrestrial vegetation canopies, coral reefs, and dense porous media, suggesting that reduced-complexity models may be successful in capturing some of the universal features of such flows 48 .
Differently from other approaches where the additional resistance introduced by canopies is modeled through a quadratic friction parametrization, we model the flow resistance generated by the canopy as a Darcy-type resistance, corrected to account for vortices penetration in the superior portion of the canopy layer through an estimate of the effective turbulent viscosity, μ t

11
. We first demonstrate that the analytical expressions of Battiato and Rubol 11 well predict the velocity profiles, volumetric flow rates and the flow resistance measured across a wide range of natural meadows with very different morphologies. The comparison focuses on profiles that do not present a secondary maximum of velocity in the lower part of the canopies, even though generalizations to this scenario are relatively straightforward and focus of current investigations. The model predictions compare well with the global resistance to flow under different conditions of plant density and morphologies.
Finally, we propose and test a novel universal scaling for the Darcian friction factor f. We demonstrate that f for canopies with low submergence is controlled primarily by the canopy permeability (1/λ) and a modified Reynolds number Re  . Importantly, we show that f/λ scales universally with the inverse of  Re , i.e. the ratio f/λ is independent of the canopy density and shape. Additionally, since for most canopies analyzed  λ = (1), we show that still exhibit a universal behavior. This result has important implications in that it allows one to quite accurately estimate the magnitude of the friction factor through measurements of discharge and slope and by means of remote sensing and acoustic techniques (to determine water level and canopy submergence). This result could pave the way to creating large scale maps of riverine vegetation-mediated friction factor estimates.
The main conclusions of this study can be summarized as follows: • Modeling submerged vegetation as a porous medium allows one to describe the flow velocity, discharge and friction factors in channels with canopy with complex morphology including grasses, woody vegetation bushes as well as prone vegetation. • Closed-form expressions derived for non-deformable porous media can be successfully adopted to describe turbulent flow over flexible canopies when (i) deflected canopy height data are available and (ii) the velocity profile does not present a secondary maximum of velocity in the lower part of the canopies due to nonuniform plant mass distribution. • For canopies with low submergence, the product between the Darcy friction factor and the canopy layer dimensionless permeability scales with the inverse of the bulk Reynolds number rescaled by the fluid to the turbulent viscosity ratio. The scaling-law solely depends on canopy height, water level and channel geometry, while it is universal across vegetation morphology. This provides a valuable tool to assess habitats sustainability associated to hydro-dynamical conditions.

Methods
Dataset used for Data Validation. The velocity profile and discharge dataset from Table 1 were collected from a variety of flume experiments with the exception of Cassan's data 56 (Cassan11, Cassan12, Cassan21, Cassan22) which investigated the flow in a real vegetated channel. The experiments included plants with morphological shapes that resembled both riverbed and riparian vegetation. Plants included grasses (e.g., eelgrass, cordgrass and barley), woody vegetation (e.g., poplars) and bushes (e.g., willows) as well as macrophytes and briophytes. The sketch of the plants morphology is included in Fig. 1, while a brief description of the plant shape is included in the following paragraph 'Dataset and Plant Morphologies' . The input parameter for the simulation -(deflected) canopy height, water level and channel slope -are listed in Table 2. Given the uncertainty in the reduced von Kármán constant for vegetated flows, data were compared with the model of Battiato and Rubol 11 using both one-fitting and two-fitting parameters. In the former we fit the canopy permeability only (black lines in Figs 2 and 3) while the von Kármán constant is set to κ = 0.19 for all the model predictions; the latter considers both the permeability and the von Kármán constant as fitting parameters (gray lines in the Figs 2 and 3). Only for submergence greater than five (i.e., datasets CassanS11, CassanS12, CassanS21, CassanS22, and Shi20), the von Kármán constant was set to κ = 0.41 in the one-parameter model. The two-fitting parameter model was then used to predict the velocity profiles of Fig. 4, while the one-fitting parameter model was used to determine the flow rate values illustrated in Fig. 5. The values of the fitting parameters and the main outputs are included in Table 2. To evaluate the error between the model predictions and the experimental data, we calculate the mean relative error for both the one-fitting and the two-fitting parameter model using the following equation: where N is the total number of data points and û i is the velocity. The datasets Shi20, Shi40 and Shi60 60 presented multiple velocity measurements at the same location: we used the averaged profile obtained by averaging the velocity values of the points at the same vertical coordinate. For Fig. 6 Table 2). The irrigation channel of Cassan et al. 56 is colonized by briophytes (Fortinalis antypiretica with stems up to few centimeters long and small leaves (1-4 mm) and long macrophytes (Ranunculus fuitanis) that were bent by the flow occupying a layer of 5-10 cm above the bed. Data are collected at two different gage stations with Station 1 mainly colonized by macrophytes and Station 2 presenting a more uniform vegetation distribution. Measurements at both stations were collected in both spring and summer. The main input and output are summarized in Table 2. Fig. 6(right). The right panel of Fig. 6 contains 164 points including flexible and rigid vegetation. We refer to the Appendix B of Poggi et al. 30 for a detailed description of the mophology used in their work.

Dataset of
Data Availability. The datasets generated during, and/or analyzed in, the current study are available from the corresponding author on reasonable request.
Buckingham Π Theorem and Universal Scaling. A standard approach to find self-similar solutions to boundary value problems is to use dimensional analysis. The drag force D can be written as a function of both geometric (K, L, H), and dynamic (Û b , μ t , μ and ρ) parameters, i.e.
where Û b is the mean (bulk) velocity, K is the canopy permeability, ρ, μ and μ t are the density, viscosity and turbulent viscosity of the flow. Since the flow is turbulent, it is expected that the dynamics is mainly controlled by ρ and μ t . Therefore, Using Buckingham Π's theorem and selecting μ t , Û b and  as repeating variables, where = + H L  is a characteristic length scale of the system, we obtain D 3 2 or, equivalently, D 4 2 In (19) and (20) are the dimensionless drag and a modified Reynolds number, defined as the product between ρ μ =Re U / b  and the fluid to turbulent viscosity ratio η = μ/μ t , respectively. Recombination of the dimensionless parameters in (20) The dimensionless parameter δλ Λ = = L K : / provides a dynamic classification between thin (Λ  1), or high submergence, and thick (Λ ≤ 1), or low submergence, canopy layers 50,69 . The former exhibit boundary layer features while the latter present a characteristic inflection point in their mean velocity profiles. Importantly, it has been shown that in thick (Λ ≤ 1) and thin (Λ  1) porous layers flow and transport are controlled either by the dimensionless height of the obstruction, δ (or its submergence) or by the obstruction dimensionless permeability λ −1 , respectively 50,69 , i.e. Since the data set considered in Figs 2 and 3 is predominantly characterized by thick porous media, i.e. small submergence (except for Cassan's profiles), we employ (25) to identify a possible scaling behaviour. We postulate the existence of self-similarity of the second kind in the form To test the proposed scaling (28) under the hypothesis that most of the resistance offered by vegetation is due to friction, i.e. ∼ C f D , we plot f/λ as a function of Re  where f is the friction factor and α = 1, see Fig. 6(left). The data analysis, discussed in the Results and Discussion Section, confirms that f/λ scales with the inverse of Re  , the modified Reynolds number for thick canopies (i.e. Λ ≤ 1, black points in Fig. 6(left)). The data regression of f/λ measured under shallow submergence (Λ ≤ 1) results in the following universal expression