Predictive model of bulk drag coefficient for a nature-based structure exposed to currents

Mangrove vegetation provides natural protection against coastal hazards like flooding and erosion. In spite of their economic and societal value, mangrove forests have experienced a worldwide decline due to human activities. Bamboo structures, formed by poles driven into the soil, are being used to create a sheltered environment for mangrove restoration. The lack of design rules for the structures has led to mixed success rates in their implementation. Improving future designs requires a better understanding of how the bamboo poles affect waves and currents. Currents cause drag forces on the poles, which depend on flow acceleration through the elements (blockage), and the distance from wakes of upstream cylinders (sheltering). We developed a model that predicts the bulk drag coefficient of dense arrays of emergent cylinders in a current, including blockage, sheltering and a balance between turbulence production and dissipation. The model could reproduce measured bulk drag coefficients from the literature within a deviation of 20%. The model also showed that anisotropic structures with small spanwise spacing and large streamwise separation maximize the bulk drag coefficient, and the energy dissipation per pole. The application of the model can guide the design of future mangrove restoration efforts.

. Incoming flow velocities U ∞ accelerate to U bl between the cylinders, an effect known as blockage. Behind the first row of cylinders, velocities reduce to U w due to sheltering effects. The relative magnitude of these effects depends on the streamwise and lateral or spanwise spacing of the cylinders ( s x and s y , respectively). www.nature.com/scientificreports/ We consequently present a physics-based model to predict the drag forces acting on emergent cylinder arrays exposed to currents, which provides a direct relationship between cylinder arrangement and c D,b . The velocities inside the arrays are estimated using a blockage factor, based on mass conservation, and a sheltering factor,  Table 1. The measurements with volumetric porosities between n = 0.92 and 0.99 were obtained from Tinoco and Cowen 22 , and the fit lines for n = 0.65-0.80 from Tanino and Nepf 21 . (b) Definition of the volumetric porosity n, given as the ratio between the fluid volume, V F = hA F over the total volume V = hA . (c) Definition of the blockage factor f b , given as the ratio between the total cross-sectional area A of the array and the constrained flow section, A c . www.nature.com/scientificreports/ based on the wake flow model developed by Eames et al. 26 . Since the model of Eames et al. 26 requires knowledge of the ambient turbulence intensity, we expand the turbulence model of Nepf 25 , including a turbulence production term by flow expansion, as done by Mossa et al. 28 , and the effects of blockage and sheltering in the wake production term. This model focuses on the local physical processes inside the structures, and it computes the bulk hydrodynamic forcing using the incoming flow velocity and flow depth as input parameters. In order to calculate the effects of the structures on the surrounding flow field (such as backwater effects or changes in the flow direction in coastal regions), the equations of the model could be built in standard free surface flow models that solve for those processes. The development of the bulk drag model is presented in the next section. Following its derivation, the model is tested against force measurements from random cylinder arrays by Tanino and Nepf 21 and Tinoco and Cowen 22 , and from regular cylinder arrays by Jansen 29 . The experiments of Jansen 29 are described in "Methods" section. The model behaviour is also explored for different cylinder configurations. Finally, the model sensitivity to different input parameters is investigated, and the model is applied to optimize future structure designs.

Model development
The analytical model consists of (1) an adapted drag formulation for closely-packed cylinder arrays, including blockage and sheltering, and (2) a turbulent kinetic energy balance, necessary to quantify sheltering. The turbulence model builds on the formulation suggested by Nepf 25 for vegetation canopies, and incorporates a turbulence production term by flow expansion, and an extended drag formulation in the wake production term. The steps to derive the equations are presented below.
Drag model. The drag forces experienced by an array of cylinders, per unit mass, can be calculated as: where c D is the drag coefficient of a single cylinder, which can be estimated using the empirical expression of White 30 , given by: where Re is the Reynolds number based on the cylinder diameter and the depth-averaged local flow velocity U. a is the projected plant area per unit volume, defined by Nepf 25 as: with d being the cylinder diameter, s the spacing between cylinders, and h the water depth. The main unknown in Eq. (1) is the local flow velocity U. If a cylinder array is sufficiently sparse, the local flow velocity could be assumed equal to the depth-averaged incoming flow velocity, U ∞ , either measured or calculated with a free surface flow model. For denser configurations, the velocity will change as the flow propagates through the array due to (1) flow acceleration between the elements (blockage), and (2) flow deceleration due to the sheltering effects of upstream rows of cylinders. Both effects are illustrated in Fig. 1c. The changes in flow velocity are included by multiplying U ∞ by a blockage factor, f b , and a sheltering factor, f s : Inserting both factors in the expression for the drag force results in Eq. (5): where the changes in velocity have been incorporated in the bulk drag coefficient, This expression provides a direct relationship between the drag coefficient of a single cylinder, c D , and bulk drag coefficients c D,b measured for cylinder arrays in laboratory experiments. Predicting the drag force thus depends on determining the values of f b and f s .
The blockage factor f b can be estimated based on mass conservation through a row of cylinders 11 , considering that the velocity will increase as the same flow discharge travels through the smaller section between the elements: where the total frontal area is A = hs y , and s y is the distance between cylinders perpendicular to the flow, centerto-center (see Fig. 1). Subtracting the frontal area of the cylinders from the total area gives the available flow area, A c : Here we are assuming that the water depth is the same just upstream and in between the cylinders. Solving for f b in Eq. (6) results in Eq. (8), see also Etminan et al. 11 : The sheltering factor f s can be estimated from the wake flow model by Eames et al. 26 , which predicts the velocity deficit behind a cylinder as a function of the distance downstream of the cylinder, s x , the cylinder diameter, the local turbulent intensity I t , and the drag coefficient: where U w is the velocity in the cylinder wake, U ∞ is the incoming flow velocity, and I t is the meant turbulent intensity, defined as I t = √ k/U ∞ 21,25 . k represents the turbulent kinetic energy per unit mass, with k = 1/2(u ′2 + v ′2 + w ′2 ) , where u ′ , v ′ , and w ′ are the instantaneous velocity fluctuations in the streamwise, lateral, and vertical direction respectively, and where the overbar denotes time averaging. The turbulent velocity fluctuations are defined as the difference between the instantaneous velocities and their mean value over a measurement period. Here we consider the depth-averaged value of the turbulent intensity, in view of the uniformity of the turbulent properties over the vertical observed inside emergent arrays 25 .
Equation (9) was developed assuming turbulent flow. Viscous effects decrease the velocity deficit 26 , with the reduction factor being given by: where Re t is the lowest Reynolds number corresponding to fully turbulent wake flow. Laminar effects are included in the wake flow model by multiplying the velocity deficit of Eq. (9) by the reduction factor f Re for Re < Re t , where the the turbulent Reynolds number is assumed equal to Re t = 1, 000 . This value is based on the observation that although a wake starts becoming turbulent at Re t ∼ 200 , drag coefficient measurements usually become constant at Reynolds numbers beyond Re t ∼ 1000 , e.g. as shown in Figure 2.7 of Sumer and Fredsoe 13 . The influence of varying Re t on the model results is investigated in "Results and discussion" section.
Defining the sheltering factor as f s = U w U ∞ , and including f Re and the bulk drag coefficient in the definition of the velocity deficit results in Eq. (11): Equation (9) also assumes that the downstream cylinder is placed inside the ballistic spreading region of the wake. The ballistic regime occurs for a downstream distance s x < L/It , where L is the integral length-scale of turbulence, and it is characterized by a rapidly decaying velocity deficit, and by a linear increase of the wake width with downstream distance. Inside the cylinder arrays, the length scale development is limited by the downstream spacing, resulting in L < s x . Considering that turbulent intensity measurements of Jansen 29 varied between I t = 0 and 0.8 inside cylinder arrays with n = 0.64-0.9, this would result in L < s x /It . This is a reasonable general assumption for the bamboo structures, since their porosity varies in a similar range. If the poles were sparsely placed, there would be a transition from ballistic to diffusive spreading of the wake. Eames et al. 26 also developed expressions for turbulent flow under the diffusive regime, which could be used in place of Eq. (9).
In the opposite case of very high pole densities, there may be a point where the elements are so closely-packed that vortex shedding is inhibited by the presence of the neighboring cylinders. Considering an analogy with a cylinder placed close to a solid boundary, vortex shedding would not take place for spanwise spacings smaller than s y /d < 1.3 13 , causing a decrease of the drag coefficient that would not be reproduced by the expression of White 30 . The application of the present model is thus restricted to s y /d > 1.3.
Balance of turbulent kinetic energy. Application of Eq. (11) requires predicting the turbulent kinetic energy. This is calculated by expanding the model developed by Nepf 25 , based on a balance between turbulence production and dissipation: where P w is the turbulent production rate and ǫ is the dissipation rate. For a dense cylinder array, k is produced by (1) generation in the wakes of the cylinders 25 , and (2) shear production by the jets formed between the elements 28 . The total turbulence production term, P w , consequently has two parts: We assume that for dense cylinder arrays these two terms are much higher than turbulence production by shear at the bed, based on observations by Nepf 25 for sparse arrays. This assumption is further tested in "Results and discussion" section.
The first term in Eq. (13), P w1 , represents turbulence production at the wakes, and can be estimated as the work done by the drag force times the local flow velocity: The second term, P w2 , represents turbulence generation due to flow expansion 28 , and can be estimated from the Reynolds shear stresses: where the energy loss due to flow expansion, E c , is modelled using the Carnot losses. Assuming that the mean kinetic energy is transformed into turbulent kinetic energy E t , and assuming isotropic turbulence, gives Eq. (17): Equation (17) enables expressing the normal Reynolds stress as a function of the incoming flow velocities and the blockage factor: The Reynolds shear stress is estimated as u ′ v ′ = Ru ′ u ′ , where the correlation factor R was given a constant value of 0.4 based on observations of Nezu and Nakagawa 31 . This value was derived for open channel flow conditions and is assumed acceptable as a first approximation, but it could vary inside a cylinder array. This is explored further in "Results and discussion" section. The velocity gradient is estimated from the velocity difference between the side of the cylinders (dominated by blockage) and the wake of the cylinders (dominated by sheltering) resulting in Eq. (19): The dissipation term, ǫ , is estimated as: The characteristic turbulent length scale l is limited by the surface-to-surface separation between the elements in the flow direction, l = min(|s x − d|, d) . This differs from the expression developed by Nepf 25 , who used the diameter as representative of the size of the eddies. We assume that in closely-packed cylinder arrays the spacing between cylinders may be smaller than the diameter, |s x − d| < d , which would limit turbulence development. The maximum value of l is set equal to the cylinder diameter. Here we also assume that for the dense cylinder arrangements, the spacing between cylinders is considerably smaller than the water depth, hence turbulence generated by bed friction is negligible.
Balancing the production and dissipation of turbulent kinetic energy results in Eq. (22): Taking the cubic root at both sides and introducing the scale factor α 1 gives Eq. (23): , which is given a default value of α 1 = 1 . The sensitivity of the model to different α 1 and R values is explored in "Results and discussion" section. k can be calculated by solving Eq. (23) iteratively, using the incoming upstream velocity U ∞ and the geometric characteristics of the structure, s y , s x , d and a, as an input. This enables determining the sheltering factor, f s = U w /U ∞ from Eq. (11). The blockage factor f b = (1 − d/s y ) −1 can also be calculated from the geometric properties of each configuration. Both coefficients can be then combined to predict the bulk drag coefficient, Deriving c D , b with the present approach relies on the assumption that the changes in water depth through the structure are small. This is a reasonable assumption given the short length of the bamboo structures in the streamwise direction, which varies between 0.7 and 1.5 m (see Fig. 1b). Longer structures that experience non-negligible changes in flow depth and velocity should be discretized, and the bulk drag coefficient should be calculated separately for the different sections. The model assumptions are discussed further in the following section.

Results and discussion
In this section we firstly present the model validation, and investigate how turbulence production and sheltering vary under different configurations. We then explore the model sensitivity to several input parameters, and finally apply the model to investigate structure design optimization.

Model validation.
The performance of the model is tested against drag measurements for regular, staggered and random emergent cylinder arrangements from the literature. A summary of the conditions tested in the different studies is shown in Table 1. The regular configurations, also denoted as in-line arrangements, consist of rows of cylinders where the downstream elements are always in one line in the streamwise direction (see configurations 1-6 tested by Jansen 29 in Fig. 7 b of "Methods" section). In the staggered arrangements, for every row the downstream elements are shifted laterally so that they are located at the center line of upstream elements, as also shown in configuration 7 of Fig. 7b. The random arrangements were obtained by distributing the cylinders using a random number generator, see Tanino and Nepf 21 . In Fig. 3 we compare the model predictions for the cases of Table 1 with two other approaches used in the literature to define the bulk drag coefficient. Figure 3a shows the bulk drag coefficient calculated from the pore velocities (based on mass conservation over the fluid volume). Figure 3b shows the drag values derived from  www.nature.com/scientificreports/ blockage factor (based on mass conservation over a cross-section, from Eq. 7). Figure 3c shows the results of the present model, which includes both blockage and sheltering effects. Using the pore velocities to estimate the bulk drag results in a general under-estimation of the drag coefficients (Fig. 3a). The blockage factor provides better estimates of the bulk drag for random arrays, but it cannot reproduce the sheltering effects observed at regular arrangements with different streamwise separations (Fig. 3b). The present model, including both sheltering and blockage, successfully reproduces the bulk drag for regular configurations, and it also provides a slight improvement for the random arrangements (Fig. 3c). The model displays a general tendency to overestimate the bulk drag of staggered and random configurations, which could be due to changes in the flow direction through such configurations.
Random and staggered arrangements have been associated to similar bulk drag coefficients in the literature 32 , which were higher than for regular configurations [33][34][35][36] . Schoneboom et al. 36 attributed the larger drag for staggered arrays to the more tortuous water flow through them. The present model assumes that the flow propagates only in the streamwise direction, and that it does not experience changes in direction. This assumption still yielded good results with the validation, especially for the densest configurations. This is expected because as the element density increases most of the total volume is occupied by cylinders. Less room for varying the spatial arrangement results in similar drag forces for regular and random arrays.
Although the model does not include changes in water level through the structures, it could still reproduce the measurements of Tanino and Nepf 21 and Tinoco and Cowen 22 , conducted with array lengths of 0.99 m and 2.84 m, respectively. This assumption may not hold for longer cylinder arrays over a fixed horizontal bed. Under those conditions the water depth could experience significant changes through the structure, which should be taken into account in bulk drag predictions. However, since the bamboo structures have a short length in the streamwise direction, such cases are beyond the scope of the present work.

Influence of spacing on hydrodynamic parameters.
Once validated, the model is applied to investigate the influence of the distance between elements on turbulence production and sheltering, and to evaluate how the previous effects translate into different c D,b values. Figure 4 shows the turbulent kinetic energy, sheltering factor, and bulk drag coefficient calculated for three values of spanwise spacing, s y /d = 1.5, 3 and 10, for streamwise separations between s x /d = [1,100].
The turbulent kinetic energy, shown in Fig. 4a, is expressed as a ratio to the turbulent kinetic energy produced by bottom friction, k o . Turbulence generation at the bed is based on the friction velocity with k o = c f ,b U 2 ∞ , where c f ,b = 0.001 corresponding to a smooth bottom. Overall, the levels of turbulence inside cylinder arrays are considerably higher than for a bare bed. The turbulent kinetic energy increases with smaller spanwise spacing s y /d , since blockage increases the drag forces, their work, and the shear production term. The largest spanwise spacing, s y /d = 10 , produces relatively lower values of k, but these are still between k = [4 − 20]k o . The turbulence levels also vary as a function of the streamwise spacing, decreasing their values for the lowest s x /d , since sheltering effects cause a strong reduction of the turbulence production terms. Higher streamwise separations reduce sheltering effects, and increase turbulence production up to a relative maximum around s x /d ∼ 2 . Beyond the maximum, the larger streamwise separations are associated to a lower number of cylinders per unit volume, a smaller projected area a, and less production of turbulent kinetic energy per unit mass.
These trends are also visible in the sheltering factor, shown in Fig. 4b, as the velocity deficit over a cylinder is inversely proportional to the level of ambient turbulence. The velocity deficit is consequently smaller for low s y /d values. Since the velocity reduction is also inversely proportional to s x /d , sheltering effects are less pronounced for higher s x /d values. This results in the bulk drag coefficients, shown in Fig. 4c, being governed by the blockage factor for s x /d > 15 , and by both sheltering and blockage for lower s x /d values. and Re t . The model sensitivity to changes around their default values is explored in Fig. 5. The scale factor α 1 is varied between 0.5 and 1.5 in Fig. 5a. The lower limit of α 1 = 0.5 is associated to a relatively low turbulence production, which in turn increases the velocity deficit on downstream elements. This results in considerable sheltering effects up to s x /d ∼ 40 . α 1 = 1 increases turbulence production and reduces the velocity deficit, causing appreciable sheltering effects up to s x /d ∼ 20 . The higher value of α 1 = 1.5 provides comparable results to α 1 = 1 . A more precise assessment of α 1 would require measurements of turbulence production and dissipation inside different cylinder configurations. Since laboratory measurements presented in the literature show that sheltering effects can be evident at a downstream distance of s x /d = 15 18 , it is concluded that α 1 = 1 provides reasonable predictions of the sheltering effect. As shown in Fig. 5b, the model results display low sensitivity to variations of the factor R, since the shear production term P w2 has a relatively lower weight on the total turbulence production in comparison with the wake production term P w1 . The influence of Re t on the bulk drag predictions is illustrated in Fig. 5c. Lower values of Re t result in stronger sheltering effects and smaller bulk drag coefficients. The largest difference between the three tested values was observed for Re = 200 , where c D,b = 2, 2.5, and 2.6 for Re t = 200, 1000, and 2000, respectively. An accurate evaluation of this threshold would require force and velocity measurements inside cylinder arrays with Reynolds numbers varying in the previous Re t range. Considering the large diameter of the bamboo poles, the Reynolds numbers in the field are most likely to be of the order of Re ∼ 10, 000 . This implies that the Re t threshold will not affect significantly the drag force predictions for the structures.

Drag maximization
The choice of pole configuration, in terms of element spacing s x /d and s y /d , is thus essential to assess the bulk drag and the resistance provided by a structure. This is conceptualized in Fig. 6. Figure 6 illustrates the computed bulk drag coefficient for different combinations of the dimensionless spacing s x /d and s y /d . The lowest value of s y /d is limited to 1.3 since, as previously discussed, below that value the expression of White 30 may not be valid. We also include solid black lines showing configurations with the same volumetric porosity. Figure 6 shows that a structure with a porosity of 80% can have an average bulk drag coefficient between c D,b = 1 and 10 depending on the element placement. The highest bulk drag coefficients are associated to rows of cylinders with a small spanwise spacing s y /d , which enhances blockage, and a large streamwise spacing s x /d , so that downstream rows experience less velocity reduction. For instance a regular structure with 80% porosity and a spanwise spacing of s y /d = 1.4 , would have a streamwise spacing of s x /d = 2.8 . This would result in a bulk drag coefficient of c D,b = 8 . If the same number of elements were placed in a uniform setting, with s x /d = s y /d = 2 , this would led to a much lower bulk drag coefficient of c D,b = 3.
Placing the rows in a staggered manner could reduce sheltering effects, but even assuming negligible sheltering, a spanwise spacing of s y /d = 2 would lead to a bulk drag coefficient of . A similar effect could be achieved with a random configuration, but predicting the net effect of the spatial changes in density on the drag would require more detailed knowledge of the cylinder density distribution. In a random arrangement the flow will tend to deflect to areas of low element density, but its trajectory will also depend on the length of the paths. A shorter path where the cylinder are more densely placed could lead to lower resistance than a longer and sparser alternative 37 . However, as previously discussed, for relatively denser structures, uniform and random arrangements should yield comparable forces. The present drag model may be implemented in large scale hydrodynamic models to evaluate the impact of currents, and the associated forces, on the cylinders. This approach would enable varying the cylinder arrangement, structure length and location, and help identify parameter combinations that optimize future structure designs. Moreover, although the present model was developed for currents, it is applicable for long waves (with KC > 100 , where KC represents the ratio of wave excursion to pole diameter) where non-stationary effects are negligible. For shorter waves (with KC < 100 ), the hydrodynamic forces also depend on additional aspects, such www.nature.com/scientificreports/ as inertial effects 38 , or turbulence enhancement by waves 39 . The influence of the previous aspects on the bulk drag coefficient is investigated further in Gijón Mancheño et al. 40 .

Conclusions
In this study a model is developed to determine the bulk drag coefficient of dense arrays of emergent cylinders, accounting for both blockage and sheltering effects. Flow acceleration through the elements (blockage) is modelled based on mass conservation through a cross-section of the array. The velocity reduction by the wakes of upstream elements (sheltering) is modelled based on the wake flow model of Eames et al. 26 , in combination with an equation to predict turbulence production by the cylinders derived by expanding the model of Nepf 25 . Turbulence production is a function of the spacing between the elements both in the spanwise and streamwise directions. Smaller spanwise spacings increase blockage and turbulence production, while smaller streamwise spacings have the opposite effect; they result in more velocity reduction on downstream cylinders, and less turbulence generation. The differences in turbulent kinetic energy also affect wake recovery, and the velocity deficit of different configurations. Higher levels of ambient turbulence result in smaller velocity deficit behind the cylinders, and less sheltering of downstream elements. When we combine blockage and sheltering effects to predict the bulk drag coefficient, we reproduce measurements from the literature for volumetric porosities between 0.64 and 0.99 within a deviation of 20%. The model also shows that reducing the lateral spacing between elements, and increasing their streamwise separation, increases the bulk drag and the dissipation per element. The application of the present model, and its development for wave flows, could help optimizing future structure designs, minimizing their material costs and erosion problems. The model may thus constitute a practical tool to increase the success of future mangrove restoration schemes.

Methods
Validation data. The data that support the findings of this study were directly obtained from the graphs of Tanino and Nepf 21 Table 1.
The locations of the instruments used during the experiments are presented in Fig. 7a. All the instruments were measuring continuously with a frequency of 100 Hz. An electromagnetic flow meter (EMF) was placed at a distance of 0.4 m upstream from the structure, at a fixed height of 0.4 m from the bottom. The EMF measured with an accuracy of 1% 41 . The instantaneous flow velocities were measured with a Nortek Vectrino acoustic velocimeter (ADV) at a fixed height of 0.4 m from the bottom. The ADV probe was installed 0.04 m upstream from the gap between two elements. The ADV measured the approaching flow before it was accelerated between two elements, and it had an accuracy of approximately 1% 42 . The output of both velocity sensors was in volts, and the velocities were obtained from linear regression, using separate calibration factors for each instrument. The hydrodynamic loads acting on one single cylinder were recorded with a SCAIME load cell mounted on the upper part of the element, measuring in volts with 0.017% accuracy 43 . The load cells were calibrated using known weights, and fitting a linear relationship between weight and voltage output. The forces were calculated by multiplying the sensor output by the calibration factor, and by the acceleration of gravity.
The bulk drag coefficients were determined by using Eq. (1) with the mean force measured at the center of each configuration, and the mean incoming velocity recorded by the EMF. The average forces and velocities were calculated using a moving average over intervals of 20 s. The top view of the structure is marked by a dashed black line, and it is illustrated for each of the tested arrangements: (C1) single cylinder with d = 0.04 m, (C2) single row with spanwise spacing between the elements of s y = 3d , (C3) single row with spanwise spacing between the elements of s y = 1.5d , (C4) multiple rows with s y = 1.5d and s x = 3d in uniform arrangement, (C5) multiple rows with s y = 1.5d and s x = 1.5d in uniform arrangement, (C6) multiple rows with s y = 3d and s x = 3d in uniform arrangement, and (C7) multiple rows with s y = 3d and s x = 3d in staggered arrangement.