Field measurement and prediction of drag in a planted Rhizophora mangrove forest

,


Introduction
The transect for drag measurement was established on September 9, 2018, a day before 169 the drag measurement was carried out. A reference tree was first identified and located at the 170 center of the transect (Fig. 1c). A visual confirmation was then made such that the above-171 ground root structures of the reference tree do not deviate largely from those of the surrounding 172 trees in terms of complexity. A 30-m long transect along the major flow direction (A-B; Fig.  173 1c), which was determined visually, was set during the ebb tide. Afterwards, the x-coordinate 174 was defined to align with the direction of mean flow during ebb tide (the major axis is shown 175 in Fig. 2d), the y-coordinate oriented laterally, and the z-coordinate oriented vertically with z = 176 0 m at the bed. The x-coordinates of A and B were defined as x1 and x2, respectively. 177 Philippines; satellite images (Google Earth) of (b) the overview of Bakhawan Ecopark, (c) 180 locations of transect A-B and the reference tree for drag measurement; and (d) photo of the 181 drag measurement site taken near the reference tree. Shoreline data are from the Global Self-182 consistent, Hierarchical, High-resolution Shorelines (GSHHG) database. 183

Measurement of vegetation variables 184
To obtain the value of the spatially averaged vegetation projected area density, a (m 2 185 m -3 ), the morphological structures of above-ground roots and stems around the reference tree 186 were extensively measured from September 13-18, 2018. Data on some trees, including the 187 reference tree, are shown in Yoshikai et al. (2021a). Ten additional trees were compiled and 188 added for a total of 23 trees for this study. The information on the locations of the measured 189 trees are provided in Fig. 2a. Some information on the vegetation parameters is provided in 190 Table 1. Here, the R. apiculata trees have multiple stems, where one tree has 3-7 stems. When 191 the main stem of a tree could not be identified in the field, the diameter of the largest stem was 192 used as DBH (diameter at breast height) of the tree. 193 m -1 10 -3.59 Scaling parameter (α 1 in Eq. (9)) --2.04 Scaling parameter (β in Eq. (9)) m 0.08 Scaling parameter (β 1 in Eq. (9)) -15.38 Root angle of approximated root shape (θ in Eq. (11)) degree -34.5 Ground level at A relative to near the reference tree m 0.045 Ground level at B relative to near the reference tree m -0.049 As described in Yoshikai et al. (2021a, b), four parameters of root were measured; these 195 are height (HR, m), horizontal distance (L, m), angle (θ, degree), and diameter (Droot, m). Then, 196 following Ohira et al. (2013), the shape of each root projected from the mean flow direction 197 was estimated from quadratic curve approximation as 198 where ψ is the azimuth angle of root to the mean flow direction; here, y = 0 at the location 200 where a root emerges, i.e., the position of y-axis varies for each root. Because the azimuth root 201 angles were not measured, a random number for ψ was given for each root in the range 0° ≤ ψ 202 < 90° using the random number generator in MATLAB. The projected area of one root can be 203 estimated by multiplying Droot with the root length provided by Eq. (2). Similarly, root volume 204 can be estimated from these parameters. By summing the projected areas of all roots in a tree 205 per vertical height, dz (m), the whole-tree root projected area per dz, aroot,i(z) (m tree -1 , here and 206 hereafter; i denotes tree index), is obtained; a value 0.01 m was used for dz throughout this 207 manuscript. 208 The stem diameter at 1.3 m height was also measured for the 23 trees. When the root 209 in a tree exceeds the 1.3 m height, the diameters above the highest root were measured. Some 210 stems branched from other stems, and in such cases, the height of the branching point was also 211 measured. Then, the whole-tree stem projected area per dz, astem,i(z) (m tree -1 ), was estimated 212 by approximating it as a patch of vertical cylinders whose stem density varies with height 213 depending on the branching of the stems. 214 Three-dimensional point clouds of the measurement site were obtained using a hand-215 site visualization (Fig. 2a). From the point clouds, locations of trees were identified, and tree 217 density of the site, ntree (tree m -2 ), was computed from the visualized tree locations (Table 1). 218 From the derived parameters aroot,i, astem,i, and ntree, the vertical profile of the parameter 219 a was calculated as 220 = 4 5677 ∑ &9 6::5,; < 9 =57>,; < , where Ntree is the number of measured trees (Table 1). Here, due to the randomness in the 222 azimuth angle in Eq. (2), the estimated value of a has some uncertainties. In this regard, the 223 value of a was computed repeatedly for 20 times, and the median value was taken as the 224 representative value of a in the area. 225

Measurement of hydrodynamic variables 226
The measurement of hydrodynamic variables for drag quantification was conducted on 227 September 10 and 11, 2018, which were spring tide conditions. A pressure sensor (U20L-04, 228 Onset Computer Corporation, USA), four electromagnetic velocity meters (EM; Infinity-EM, 229 JFE Advantech, Japan), and one Acoustic Doppler Velocimeter (ADV; Nortek Vector, 230 Norway) were deployed around the reference tree ( Fig. 2a) for the two-days measurement. The 231 EMs were deployed to measure the near-bed velocities behind (P1), in front (P3), and the sides 232 (P2, P4) of the reference tree relative to the flow direction during the ebb tide (Fig. 2a, b). The 233 body of EM was buried in the mud to position the probe at 5 cm above the bed as done in 234 Schettini et al. (2020). The ADV was deployed around 3 m away from the reference tree in a 235 downward-looking orientation (Fig. 2a, c), where the center of the sampling volume was placed 236 at 5 cm above the bed for bed drag quantification (Pope et al., 2006). The ADV was set to 237 collect data with 16 Hz for 1 minute (960 samples) every 10 minutes. The pressure sensor and EMs were also set to collect data every 10 minutes. Using data from the deployed EMs, it was 239 confirmed that the flows had a distinct axis and did not rotate during the tidal cycles (Fig. 2d).  The water level differences between the transect ends A and B, Δη (m), were measured 250 during the ebb tides using the water leveling method. A schematic of the measurement setup is 251 shown in Fig. 3a. The method is based on the principle that the water level at both ends of a 252 conduit will equalize based on atmospheric pressure. Plastic poles were installed at the transect 253 ends A and B and a 35-m long plastic tube (inner and outer diameters: 6 and 8 mm, 254 respectively) spanned between the poles as illustrated in Fig. 3a. Also, a 1.5-m long plastic 255 tube with the same inner and outer diameters was placed onto pole B vertically (Fig. 3a, b). 256 When the ground was submerged during flood tide on a measurement day, the water 257 connectivity within the tubes was ascertained by removing any trapped air from the upward-258 oriented tube end using a syringe. This made the water level inside the tube equalized at a 259 location where the downward-oriented tube end is placed, hence the water levels at A (η(x1)) 260 and B (η(x2)) were made visible at B (Fig. 3c). The downward orientation of the tube end was 261 to prevent the effects of pressure created by flows on the water level inside the tube. When the 262 water became still during high tides, the same level of η(x1) and η(x2) was confirmed. During 263 ebb tides, Δη was measured using a caliper with 0.1 mm resolution every 10 minutes 264 synchronized to the timing of data collection by the deployed sensors (Fig. 3c). The Δη was 265 recorded as 0 mm when the water level difference was too small to measure even though the 266 waters were flowing. 267 The water depths at A, B, and near the reference tree were measured manually when 268 the water was still at high tide. From these water depths, the ground levels at A and B relative 269 to the site near the reference tree were calculated (Table 1). 270 In conjunction with the sensor data collection and the water level difference 271 measurement, vertical profiling of flow velocity was conducted at the four locations around the 272 reference tree (P1-P4; Fig. 2). An electro-magnetic current meter equipped with a pressure 273 sensor (AEM213-DA, JFE Advantech) was used for the profiling. The sensor is connected to 274 a display unit with a cable and collects data at 1 Hz; one person stood on the root system of the 275 reference tree above the water surface and slowly moved the sensor down (~1.0 cm s -1 ) from the water surface to the bottom using a cable. When the water became shallower than around 277 20 cm, a propeller velocimeter (CR-11, Cosmo Riken, Japan) was used for the profiling instead 278 of the AEM213-DA. Its small propeller (~ 2 cm diameter) made it possible to measure velocity 279 within a thin layer and is well-suited to profile shallow waters. In this case, the flow velocities 280 along the transect (x-axis) at the surface, middle, and bottom layers were measured at the four 281 locations around the reference tree.

Data processing and bed drag estimation 289
To obtain the streamwise mean flow velocity profile, the velocities along the x-axis, 290 which were measured by the AEM213-DA, were binned using 0.05-m depth-width and averaged in each bin. The channel mean velocity, U, was calculated by averaging the mean 292 velocities in the bins (or layers, in the case of data from the propeller velocimeter) of the four 293 locations (P1-P4), based on the assumption that the average of the four locations represents 294 the spatial average in the area. To check the validity of this assumption, a particle tracking 295 velocimetry (PTV) survey was conducted around the reference tree in March 2019. See Text 296 S1 and Fig. S1 for the details. 297 The bed drag, Fbed (m 2 s -2 ), was quantified from the measured Reynolds stress provided 298 by the ADV data (see Text S2 for the details). A bed drag coefficient (Cbed) was then 299 determined by fitting the measured Fbed and U in following equation of the quadratic drag law 300 where the value of Cbed was determined as 4.2 × 10 -3 with R 2 = 0.55 (see

Estimation of vegetation drag and drag coefficient 305
Drag by vegetation was quantified from the depth-averaged momentum balance. The 306 inertial terms were significantly small compared to the pressure gradient (more than 20 times 307 smaller), and thus they were neglected as done in other works (Nepf, 1999 where g is the gravitational acceleration (m s -2 ), h is the water depth (m), and Fveg,obs is the 312 vegetation drag (m 2 s -2 ) quantified from field data. Here, we assumed that the water flux is 313 where Uref and href are channel mean velocity and water depth near the reference tree, respectively, and the bed slopes between 315 A-reference tree and B-reference tree are constant. Then, following Lentz et al. (2017), Eq. 316 (5) can be rearranged by horizontally integrating between x1 and x2 as 317 where the angle bracket denotes the spatial average between the transect ends A-B and 319 This equation gives estimates of the mean vegetation 320 drag between the transect, 〈 , CN 〉. Similarly, by assuming that value of CD does not vary 321 between the ends of the transect, integration of the drag model Eq. (1) between x1 and x2 yields 322 The value of CD was derived by equating Eqs. (6) and (7) for each measurement.
where HRk and HRmax are the root heights (m) of k th highest and the highest root in a root system, 332 respectively, and S is a scaling factor. The parameters S and HRmax can be expressed as 333 functions of DBH, and thus, HRk is a function of DBH as 334 where β, β1 and α, α1 are the scaling parameters for HRmax and S, respectively. The values of 336 these parameters for our study site were derived in Yoshikai et al. (2021a) (Table 1). Similarly, if the t th highest root is the one with the minimum height in a root system, t is the largest integer 338 number that satisfies 339 where HRmin is the critical height (m) to be given as a model parameter. From Eqs. (9), (10), 341 the vertical variation of root density per dz in a tree is modeled. 342 To compute the root projected area from the modeled root density, an empirical 343 relationship of nroot,i(z) and aroot0,i(z) provided in Fig. S3 was used, where nroot,i(z) is the number 344 of roots per dz in tree i (root m -1 tree -1 ), and aroot0,i(z) is the root projected area per dz with zero 345 azimuth angles (m tree -1 ). The strong linear relationship between nroot,i(z) and aroot0,i(z) suggests 346 that the individual roots can be approximated to a single linear shape assuming a uniform root 347 diameter as 348 If Eq. (11) is applied, the slope of the nroot,i(z) and aroot0,i(z) relationship stands for Drootdz/sin(-350 θ). By applying the average root diameter (Droot = 0.03 m; Table 1), the value of the root angle 351 θ was determined as -34.5° for our study site ( Fig. S3a; Table 1). As with the field data, random 352 numbers were given to ψ in Eq. (11) for each root. The parameter a was then calculated from 353 the modeled root projected area using Eq. (3) for 20 times, and the median value was taken as 354 the representative value of model prediction; the observed value was used for astem,i(z), which 355 can be easily measured in the field.    and 11, 2018. Note that some data on Uref and Δη are absent due to instrument problems and 386 measurement setup maintenance, therefore the number of measured 〈 , CN 〉 is smaller than 387 those of Uref or Δη (Fig. 5); also, 〈 , CN 〉 was not derived when the Δη recorded was 0 mm. 388

Results
The Uref was generally around 1.5 times larger than the velocities near the bed measured by 389 EM sensors but became comparable when the water depth decreased (href < 0.2 m). While these 390 patterns were consistent during the two-days measurement, the Uref on September 10, 14:00, 391 was smaller than the velocity from EM sensors, possibly due to an unreliable measurement of 392 Uref using the propeller velocimeter (Fig. 5b). Velocity magnitude measured by ADV during flood tide was significantly lower than EM-measured velocity or Uref (Fig. 5b, f), probably due 394 to local influence of nearby roots at the upstream side (Fig. 2c). 395 The variations of the measured Δη were 1.2-11.4 mm (Fig. 5c, g). The Δη increased as 396 water depth decreased. Accordingly, the vegetation drag per water volume 〈 , CN 〉 ℎ c d ⁄ 397 showed an increase with decreasing water depth (Fig. 5d, h). The bed drag per water volume 398 〈 C 〉 ℎ c d ⁄ was significantly small compared to the vegetation drag 〈 , CN 〉 ℎ c d ⁄ , more 399 than 15 times smaller during most of the measurement time. of the spatially averaged velocity (u) also showed significant decrease below 0.25 m (Fig. 6b), 414 corresponding to a significant increase in a (Fig. 4). The profile of u showed agreement with a 415 theoretical predictor of spatially averaged velocity profile (red line in Fig. 6b) (Fig. 7a). Note that the data taken on September 11, 14:00, was 429 excluded from the line fitting as the data of Uref may not be accurate (see Fig. 5b). Instead, the 430 〈 〉 and drag averaged for unit vegetation projected area 〈 , CN 〉 〈 〉 ⁄ showed higher 431 correlation (R 2 = 0.901), and separate line fittings did not show significant difference in the 432 line slopes (Fig. 7b).  Fig. 4. The result showed a high coefficient of 457 determination (R 2 = 0.86) for the vegetation drag averaged for unit water volume (Fig. 9).  Fig. 4. 462

Prediction of drag using the Rhizophora root model 463
The Rhizophora root model well-predicted the overall vertical profile of a composed 464 of multiple order roots, using a parameter setting of θ = -34.5°, a value determined for the 465 Bakhawan Ecopark study site (Fig. S3a), and HRmin = 0.01 m (Fig. 10a). The modeled 466 vegetation drag computed with the modeled a and CD = 1.0 showed good agreement with the 467 measured drag, with a slope of 1.06 and R 2 = 0.84 of the linear fitted line (Fig. 10b). The use 468 of θ value obtained in another mangrove forest (-41.9°; Fig. S3b) resulted in the 469 underestimation of a due to the steeper angle of the approximated root shape (Eq. (11)) 470 specifically at the lower part (z < 0.3 m) (Fig. 10c). Due to the underestimation of a, the 471 predicted vegetation drag also showed underestimation trend, with the fitted line slope of 1.18, 472 while the R 2 value did not vary significantly compared to when θ = -34.5° (Fig. 10d). 473 The increase in the value of HRmin from 0.01 m to 0.1 m resulted in underestimation of 474 the vegetation drag as seen in the increased slope of fitted line (Fig. 10e)

Flow and drag in the studied field mangrove forest 489
The spatially uniform distribution of the Rhizophora trees with the same age at the site 490 investigated in this study (Bakhawan Ecopark in Aklan, Philippines; Fig. 1d; Fig. 2a) represents 491 a setting that previous laboratory-based studies have examined using model mangroves (Zhang 492  s -1 at the maximum (Fig. 5b, f), which is also comparable to the velocity measured in other 508 mangrove forests (e.g., Chen et al., 2016;Horstman et al., 2021). The mangrove forest 509 investigated in this study is thus considered to have a typical vegetation projected area density 510 and tidal flow regime that can be observed in other mangrove forests. This implies that the 511 insights obtained in this study are applicable to other mangrove forests. 512 The normalized local velocity (up/u0.1) showed larger spatial variations at higher 513 elevation (z > 0.25 m) compared to lower elevation (Fig. 6a). Generally, roots are more 514 clumped around the stem at the higher part of the root system, making locally low root blockage 515 areas especially at the sides of a tree (P2, P4; Fig. 1d). The relatively higher velocity at P2 and 516 P4 may be due to flow redistribution to the low blockage area (Maza et al., 2017), and the 517 lower velocity at P1 or P3 may be due to the influence of wakes generated by roots and stems 518 or velocity deceleration upstream of the clumped roots (Chen et al., 2012). Roots are spread 519 widely at the lower part of the root system (Méndez-Alonzo et al., 2015; Fig. 1d) making a 520 relatively uniform root distribution, which may explain the smaller spatial variations of 521 velocity at lower elevations. 522 The profile of normalized velocity averaged for P1-P4 (u) showed a good agreement 523 with the theoretical model of Lightbody and Nepf (2006), which predicts the profile of spatially 524 averaged velocity in vegetations with vertically-varying frontal area (Fig. 6b). This model is 525 One key feature observed in the field mangrove forest is the depth dependence of drag 535 per water volume as seen in Fig. 7a. The different slopes of the relationship between 〈 〉 and 536 〈 , CN 〉 ℎ c d ⁄ depending on the water depth indicates the enhancement of drag per water 537 volume relative to flow velocity when the water depth decreased. This may be considered the 538 result of the vertical variation of the parameter a, leading to increased depth-averaged 539 vegetation projected area to exert drag per unit water volume as the water depth decreases (Fig.  540   4). This highlighted the difference of drag characteristics in Rhizophora mangrove forests from 541 an array of vertical emergent cylinders and the difficulty in parameterizing mangrove drag 542 effects in roughness parameters such as Manning's roughness coefficient and Chezy coefficient 543 with 〈 , CN 〉 〈 〉 ⁄ (Fig. 7b), suggesting that the drag exerted per unit vegetation area solely 546 depends on the square of flow velocity. This signifies that the quadratic drag law is applicable 547 to the studied field mangrove forest. indicates turbulent wake structures (> 1,000) throughout a tidal phase, and the derived CD is 553 independent of Re (Fig. 8). Interestingly, the CD derived for the studied field mangrove forest 554 also showed values around 1.0, close to the ones obtained for the model mangroves despite the 555 complicated root systems that field mangroves have. This value also agrees with the value (1.0) 556 which is typically used for the drag coefficient of other type of vegetation (e.g., seagrass) at 557 This study used spatially averaged equations (Eqs. (1), (4)-(7)) for deriving CD. 562 Therefore, the estimates of CD could be significantly biased by the error in measuring the 563 channel mean flow velocity, U. While it is challenging to obtain the true value of channel mean 564 velocity and assess the measurement error, we refer to the results of the PTV survey conducted 565 around the reference tree (Text S1, Fig. S1). The results suggest that the velocity averaged for 566 the four locations (P1-P4) deviates 10% to 20% from the PTV-estimated spatially averaged 567 velocity. This deviation leads to CD error estimates of 20% to 35%, which are close to the 568 variations of the derived CD values (Fig. 8). We thus consider that the derived CD in the field 569 mangrove forest may have errors of approximately 20-35% and the variations of the obtained 570 CD are attributed to the errors in measuring the channel mean flow velocity. 571 Our observation of the quadratic dependence of drag on velocity (Fig. 7b) and the 572 obtained value of CD ≈ 1 (Fig. 8) suggests the applicability of the drag model, Eq. (1), to field 573 mangrove forest settings. The good agreement with the modeled and observed drag shown in 574 Further field-based studies are needed to consider these aspects. 601

Implications for representing mangrove drag effects in hydrodynamic models 602
Representing mangrove drag effects using Eq. (1) in hydrodynamic models have been 603 challenging because of the need for information on vegetation morphology for the parameter 604 a, which is labor-intensive to obtain in the field. This study presented a measure to predict a in 605 addition to the field estimates of CD, which none of the previous field-based studies on 606 mangrove drag were able to consider (e.g., Mazda et al., 1997;Horstman et al., 2021). We used 607 the Rhizophora root model of Yoshikai et al. (2021a) to predict a, which is based on the 608 allometric scaling of root structures. This model is valid for complicated root systems 609 composed of multiple order roots, and accurately predicted the vertical profile of a in our study site (Fig. 10a). The good agreement of the modeled drag using the field-derived CD ≈ 1 and the 611 predicted a with the measured drag (Fig. 10b) suggests that the drag in Rhizophora mangrove 612 forests in the field can now be predicted once the input parameters of the Rhizophora root 613 model are given. Note that because the roots higher than the 1 st order could dominate in a 614 specifically at lower elevations as shown in Fig. 4 Some considerations should be noted when using the Rhizophora root model, especially 625 on its parameter settings. First, the scaling parameters of root systems (Table 1)  (results not shown). Therefore, we suggest the measurement of the minimum root heights in 634 the field to find a representative value of HRmin at the site in addition to the parameters for obtaining the scaling relationships. Lastly, the root angle of the approximated linear root shape 636 in Eq. (11) seems to vary depending on the site or species (Fig. S3). The use of root angle 637 determined for Fukido mangrove forest, which is 7.5° steeper than our study site, affected the 638 prediction of the a and the drag to some extent as shown in Fig. 10c-e. The root angle was 639 determined from the relationship between nroot,i(z) and aroot0,i(z), and both parameters are labor-640 intensive to obtain in the field. Hence, determining the representative root angle for the site of 641 model application may be challenging. Nevertheless, the responses of the predicted a and drag 642 with the different parameter settings provided in Fig. 10 can be used as benchmark for model 643 uncertainty when applying the settings to other mangrove forests. Notably, the drag can still be 644 predicted with reasonable accuracy using estimates of root angle from the other mangrove 645 forests (Fig. 10c-e), thus highlighting the significance of this work in contributing to a better 646 prediction of drag in mangrove forests. 647

648
This study presents the drag force and drag coefficient estimated from a 17-year-old 649 planted Rhizophora mangrove forest based on comprehensive hydrodynamics and vegetation 650 morphology data collected from the field. The vegetation projected area density, a, showed 651 nearly exponential increase towards the bed mainly due to root branching, highlighting the 652 complex root systems of mangroves. Consequently, the drag averaged for unit water volume 653 showed depth dependence relative to velocity, suggesting the difficulty in parameterizing the 654 drag effects of Rhizophora mangroves using bed roughness parameters. Instead, the drag 655 Ferrera for providing language help. We are also grateful to Mr. Allan Quimpo of the Kalibo 679 Save the Mangroves Association (KASMA) and the staff of Bakhawan Ecopark.

Text S2. Acoustic Doppler Velocimeter (ADV) data processing
The velocity data collected by the ADV were despiked using the phase-space method described in Mori et al. (2007). The despiked velocities (eastward, northward, and vertical) were rotated to give the velocities along the x, y, and z-axes, where the instrument tilt was corrected to make the averaged vertical velocity zero (Lee et al., 2004). Bed drag (F bed , m 2 s -2 ) was then determined from Reynolds stress, (− ′ ′ ), where u' and w' are the velocity fluctuations of xand z-axis components (m s -1 ), respectively, and the overbar denotes the time average (note that velocities in the equations in the manuscript denote timeaveraged values without the overbars).
As shown in Figs. 5b and 5f, the velocity measured by ADV during the flood tides largely deviated from the EM-measured velocity and U, possibly due to the local influence of nearby roots. Thus, the Reynolds stress measured during flood tides might have been affected by the wakes generated by the roots aside from the bottom friction. In this regard, we excluded the data during flood tides when estimating the bed drag coefficient.
The estimated drag coefficient, C bed , was 4.2 × 10 -3 (Fig. S2). This value is higher but in the same order of magnitude as the drag coefficient observed in muddy tidal environment (e.g., 2.5 × 10 -3 in Mariotti and Fagherazzi, 2012).