**Picoplankton Distribution and Activity in the Deep Waters of the Southern Adriatic Sea**

#### **Danijela Šanti´c 1,\*, Vedrana Kovaˇcevi´c 2, Manuel Bensi 2, Michele Giani 2, Ana Vrdoljak Tomaš 1, Marin Ordulj 3, Chiara Santinelli 2, Stefanija Šestanovi´c 1, Mladen Šoli´c <sup>1</sup> and Branka Grbec <sup>1</sup>**


Received: 19 July 2019; Accepted: 8 August 2019; Published: 10 August 2019

**Abstract:** Southern Adriatic (Eastern Mediterranean Sea) is a region strongly dominated by large-scale oceanographic processes and local open-ocean dense water formation. In this study, picoplankton biomass, distribution, and activity were examined during two oceanographic cruises and analyzed in relation to environmental parameters and hydrographic conditions comparing pre and post-winter phases (December 2015, April 2016). Picoplankton density with the domination of autotrophic biomasses was higher in the pre-winter phase when significant amounts of picoaoutotrophs were also found in the meso-and bathy-pelagic layers, while *Synechococcus* dominated the picoautotrophic group. Higher values of bacterial production and domination of High Nucleic Acid content bacteria (HNA bacteria) were found in deep waters, especially during the post-winter phase, suggesting that bacteria can have an active role in the deep-sea environment. Aerobic anoxygenic phototrophic bacteria accounted for a small proportion of total heterotrophic bacteria but contributed up to 4% of bacterial carbon content. Changes in the picoplankton community were mainly driven by nutrient availability, heterotrophic nanoflagellates abundance, and water mass movements and mixing. Our results suggest that autotrophic and heterotrophic members of the picoplankton community are an important carbon source in the food web in the deep-sea, as well as in the epipelagic layer. Besides, viral lysis may affect the activity of the picoplankton community and enrich the water column with dissolved organic carbon.

**Keywords:** picoplankton community; deep-sea; Southern Adriatic; Mediterranean Sea

#### **1. Introduction**

Autotrophic members of the picoplankton community are important primary producers, while bacteria are consumers of the dissolved organic matter that can originate from primary production, thus transferring carbon towards higher trophic levels. Bacteria represent the main component of plankton, which is involved in the degradation of organic matter and in transforming inorganic compounds into the shapes adequate for primary producers. The role of the picoplankton community has become more important in oligotrophic and phosphorus-limited (P-limited) areas, such as the open sea area of the Adriatic Sea (Mediterranean Sea), where bacteria, *Synechococcus*, *Prochlorococcus,* and picoeukaryotes play an important role in the production and flow of biomass and energy [1–3] than in eutrophic areas. However, previous studies on picoplankton communities were mostly focused on investigating the epipelagic layer (i.e., depths less than 200 m). The deep-sea is characterized by the absence of light, i.e., conditions that are unfavorable for the primary production. Tanaka and Rassoulzadegan [4] pointed out the importance of bacteria and their biomass in carbon flux in the deep-sea. Moreover, Arístegui et al. [5] have highlighted that the deep ocean represents a key site

for remineralization of organic matter and long-term carbon storage. The discovery of cyanobacteria *Synechococcus* in the deep part of the Adriatic Sea revealed that they could be used to gain a better understanding of the effects of deep-ocean convection, such as ventilation and renewal of deep waters [6]. Hence, the vertical distribution of the picoplankton in the open southern Adriatic Sea, below the euphotic zone, has recently started to be investigated more intensively [7–10]. Aerobic anoxygenic phototrophs (AAPs), a bacterial group recently discovered in the Adriatic Sea [11,12] and primarily unknown in the investigated area, are up to 3 × bigger than other bacterial cells [13,14], and hence they could represent a remarkable source of carbon in the marine environment. These photoheterotrophic microorganisms can harvest light energy using pigment bacteriochlorophyll *a* to supplement their primarily organotrophic metabolism only in the presence of oxygen [15].

The Adriatic Sea is an elongated semi-enclosed basin of the Eastern Mediterranean Sea. According to its morphology and bathymetry, it can be divided into three sub-basins (northern, middle, and southern). The Southern Adriatic Pit (hereafter SAP) is the deepest area, with a depth reaching ~1250 m. Adriatic is characterized by a cyclonic basin-scale circulation. Through the Strait of Otranto at its southern end (~80 km wide, with a sill depth of ~800 m) [16,17], the Adriatic exchanges water and mass with the adjacent Ionian Sea. The Eastern Adriatic Current (EAC) brings northward warm and saline waters from the Ionian Sea and the Western Adriatic Current (WAC), transporting southward colder and fresher waters [18]. Waters from the Ionian Sea enrich the P-limited Adriatic Sea [19,20] with nutrients and organic substances, causing changes in the food web [21] and the distribution of organisms [6,8,10,22]. The prevalent heat and salt import from the Ionian involve the surface and intermediate layers of the Adriatic Sea, between the surface and 800 m depth. According to Gaˇci´c et al. [23], periodical changes in the sense of the rotation between cyclonic and anticyclonic phases of the North Ionian Gyre (NIG) deviate the branch of the Atlantic water, entering the Eastern Mediterranean through the Strait of Sicily, towards the Adriatic (during anticyclonic phase) or towards the Cretan and Levantine basins (during cyclonic phase). Under these two opposite conditions, changes in the heat and salt content of the SAP may occur out of phase compared to those observed in the Levantine basin [24]. Moreover, water masses flowing into the Adriatic along its eastern flank have important ecological implications. Indeed, the pool of available nutrients in the SAP changes according to the phases of the NIG [25–27]. The deep part of the SAP occasionally receives the North Adriatic Dense Water (NAdDW), that ventilates the deepest layers after winter cooling [28,29]. Locally, the open-sea winter convection triggers vertical mixing and produces dense water (Adriatic Deep Water, AdDW) that reaches depths between 300 and 1000 m [30–32]. Subsequently, AdDW outflows across the Strait of Otranto sinks into the Ionian Sea abyss and ventilates the deep layers of the eastern Mediterranean [33].

The picoplankton community is exposed to sudden physical-chemical changes in a dynamic environment. The ability of certain members to acclimate physiologically determine their presence and activity or absence in the water column. Hence, environmental conditions (vertical convection, lateral exchanges, long-term trends, and multiannual oscillatory signals) make the Southern Adriatic an ideal laboratory to study vertical distribution and activity of picoplankton members in the deep zones (below epipelagic layer), where they can represent a significant source of carbon through their biomass.

The primary goal of the ESAW (evolution and spreading of the Southern Adriatic Waters) cruise was to assess and compare the pre-and post-winter hydrographic and environmental (physical, biogeochemical) conditions in the middle and southern Adriatic basins. The objective of the present study was to evaluate for the first time, the abundance of aerobic anoxygenic phototrophs (AAPs) in the deep waters of the Southern Adriatic Pit (hereafter SAP). Beside the AAPs, we determined the distribution of *Synechococcus*, *Prochlorococcus*, picoeukaryotes, and heterotrophic bacteria and estimated their contribution to carbon budget from epipelagic to deep waters. Our study shows and discusses results obtained from the analysis of picoplankton biomasses and their activity concerning environmental parameters and hydrographic conditions and the effect of grazing and viral lysis.

#### **2. Material and Methods**

#### *2.1. Physical and Biochemical Sampling*

In this study, we used data collected during the two ESAW cruises (Evolution and spreading of the Southern Adriatic Waters), supported by the EUROFLEETS2 program. The first one was conducted in December 2015 (10–15 December). The second one was conducted in April 2016 (5–10 April), when the hydrographic and environmental conditions were supposed to be changed after the wintertime occurrence of vertical convection, and consequent dense water formation. The intention was to measure the signals of a possible spread of dense water from the middle to the southern Adriatic, because the former may be both an accumulation site and an occasional formation site for dense waters [34]. Here, we focused only on data collected in the SAP, at hydrological stations regularly spaced along a transect between Bari (Ba; Italy) and Dubrovnik (Du; Croatia, Figure 1). A conductivity-temperature-depth (CTD) SBE 911plus probe (Sea-Bird Electronics, Bellevue, WA, USA) measured vertical profiles of temperature (T), conductivity (C), dissolved oxygen (DO), fluorescence, and turbidity. The sampling rate was set to 24 Hz. The T-C sensors were controlled and calibrated before and after the cruise at the Calibration facility center at OGS (Trieste, Italy). For water sampling purposes, the CTD probe was coupled with a rosette sampler, holding 12 Niskin bottles (8-L capacity). Salinity (S), potential temperature (θ), and potential density anomaly (σθ) were calculated using the MATLAB toolbox TEOS-10 (http://www.teos-10.org/software.htm) and Ocean Data View software [35]. The reference pressure for θ and σ<sup>θ</sup> was 0 dbar. Data were processed and quality checked according to MyOcean in situ quality control standards and methodology. The final vertical profiles consisted of the data averaged every 1 dbar, with overall accuracies ±0.005 ◦C for T, and ±0.005 for S, and 2% of DO saturation. Due to the malfunctioning of the SBE43 sensor for DO detection, these data were missing at stations 8 and 9 during the April cruise. The missing profiles were substituted by the data from the closest stations, namely, 7 and 10 (Figure 1), assuming that the main pattern of the vertical structure was similar in the deep SAP region (as compared in Figure 2). Sampling depths for the nutrients and dissolved organic compounds were 2, 45 (56 m in April), 84 (100 m in April), 220 (204 m in April), 400, 500, 662 (600 m in April), 800, 940 (900 m in April), 1000, 1100, 1140 m at station 8; 2, 50 (62 m in April), 96 (100 m in April), 193 (201 m in April), 400 (350 m in April), 520 (missing in April), 700, 800, 900, 1000, 1100, 1208 m at station 9. Dissolved organic carbon (DOC) samples were occasionally sampled at a smaller rate. During sampling, there were no pieces of evidence on neither water leakage from the Niskin bottles nor some other possible source of contamination of the water samples. Moreover, there were no striking patterns among measured parameters which could have indicated such kind of a problem. Therefore, we are confident that the samples were taken correctly.

Water samples for measuring dissolved inorganic nutrient concentrations (nitrite, NO2, nitrate, NO3, ammonium, NH4, phosphate, PO4, and silicic acid, H4SiO4) were prefiltered through glass-fiber filters (Whatman GF/F) immediately after sampling and stored at −20 ◦C in polyethylene vials until performing analysis in the on-shore laboratory. The samples were defrosted and analyzed colorimetrically with a QuAAtro Seal Analytical continuous segmented flow analyzer, according to Hansen and Koroleff [36]. Detection limits for the procedure were 0.003 μM, 0.01 μM, 0.04 μM, 0.02 μM, and 0.02 μM for NO2, NO3, NH4, PO4, and H4SiO4, respectively. The accuracy and precision of the analytical procedures at low concentrations were checked periodically through the quality assurance program QUASIMEME, and the relative coefficient of variation for five replicates was less than 5%. Internal quality control samples were used during each analysis. As in the case of inorganic nutrients, samples for dissolved organic nitrogen and phosphorus analysis were filtered through Whatman GF/F glass fiber filters, collected in acid-washed polyethylene vials rinsed with seawater, and frozen immediately (−20 ◦C) until laboratory analysis. Total dissolved inorganic nitrogen (TDN) and phosphorus (TDP) were determined after quantitative conversion to inorganic N and P by persulfate oxidation [36] and subsequent analyses of nitrate + nitrites and phosphates. Dissolved organic nitrogen (DON) and phosphorus (DOP) were computed from the relationship DON = TDN

– (N-NH4 + N-NO2 + N-NO3) and DOP = TDP − P-PO4. DOC samples were filtered immediately after collection through 0.2 μm membrane filters (Sartorius, Minisart, SM 16534, Goettingen, Germany) and stored at 4 ◦C in the dark until analysis (within 3 months). DOC samples were analyzed using a Shimadzu TOC-VCSN. Samples were acidified with 2N HCl and sparged for 3 min with CO2-free pure air, to remove inorganic carbon before high-temperature catalytic oxidation. One hundred microliters of the sample were injected into the furnace after a four-fold syringe washing. From 3 to 5, replicate injections were performed until the analytical precision was within <2% (±1 μM). The calibration curve was determined by using four different solutions of potassium phthalate, in the same concentration range as the samples. The reliability of the measurements was checked daily using consensus reference materials (CRM) kindly supplied by Prof. D. A. Hansell, University of Miami. At least two CRM and two low carbon water (LCW) analyses were performed for each analytical day (measured value = nominal value ± 0.5 μM).

**Figure 1.** The Adriatic Sea. Blue dots indicate CTD (conductivity-temperature-depth) stations of the ESAW (Evolution and spreading of the Southern Adriatic Waters) cruises. Red line indicates the SAP (Southern Adriatic Pit) transect. Biological sampling concerns stations 8 and 9.

**Figure 2.** Vertical distribution of (**a**) potential temperature, (**b**) salinity, (**c**) dissolved oxygen, along the Bari-Dubrovnik (Ba-Du) section during the ESAW1 cruise in December 2015. Grey dots indicate locations of picoplankton sampling. Panels (**d**–**f**) show the vertical profiles of the same properties at stations 8 and 9.

#### *2.2. Picoplankton Analysis*

Sampling depths for the picoplankton analysis were 2, 45 (56 m in April), 220 (204 m in April), 400, 662 (600 m in April), 800, 1000, 1140 m at station 8; 2, 50 (62 in April), 96 (100 m in April), 193 (201 m in April), 400 (350 in April), 520 (missing in April), 700, 800, 1000, 1208 m at station 9. Flow cytometry was used to determine the abundances of *Synechococcus, Prochlorococcus*, picoeukaryotes, and heterotrophic bacteria [37]. Samples for autotrophic cells analysis (2 mL) were preserved in 0.5% glutaraldehyde, frozen at −80 ◦C, and stored until analysis (5–10 days). Samples for analysis of bacteria were preserved

in 2% formaldehyde and frozen until analysis (5–10 days). Autotrophic cells were divided into two groups: cyanobacteria (*Synechococcus* and *Prochlorococcus*) and picoeukaryotes, distinguished according to light scattering, cellular chlorophyll content, and phycoerythrin-rich cells signals, respectively. Bacterial abundance was determined in scatter plots of particle side scatter versus Sybr Green I fluorescence related to cellular nucleic acid content, to discriminate bacteria from other particles. According to the cellular nucleic acid content, the bacterial population is divided into two sub-groups, HNA (High Nucleic Acid content) and LNA (Low Nucleic Acid content) bacteria. Abundances of Sybr Green-I-stained heterotrophic nanoflagellates (HNF) were also determined by cytometry [38]. For obtaining abundances, samples were processed on a Beckman Coulter EPICS XL-MCL with a high flow rate from 1 to 1.2 μL s−1. AAPs were sampled only in April and were determined using the protocol described by Mašìn et al. [39]. Cells were collected on 0.2-μm polycarbonate (PC) filters by filtration and dyed with 4 ,6-diamidino-2-phenylindole (DAPI) using 3:1 mixture of Citifluor™ AF1 and Vectashield® after drying. AAP bacteria were enumerated using an Olympus BX51 microscope equipped with an Olympus UPlanSApo 100×/1.40 OIL, IR objective, and software for image analysis (CellSens, Münster, Germany). The microscope was equipped with a Hg Lamp U-RFL-T, Olympus, for excitation. Fluorescent images were taken using an XM10-IR camera. Three epifluorescent filter sets were used, DAPI, IR, and chlorophyll, to create the composite image. These images were afterward used for distinguishing between organisms that contain bacteriochlorophyll *a* and chlorophyll *a*, but also for determining the number of heterotrophic bacteria, cyanobacteria, and AAP bacteria in each sample. The abundance of virus-like particles (VLP) was determined, as described in Noble and Fuhrman [40]. Collected samples were preserved in formaldehyde (2%, final concentration), flash-frozen in liquid nitrogen, and stored at −80 ◦C until analysis, which was performed in the laboratory immediately after the end of the cruise. Preserved samples (2 mL) were filtered through 0.02-μm pore-size filters (Anodisc; diameter: 25 mm; Al2O3, Whatman, Maidstone, UK) and stained with Sybr Green I (stock solution diluted 1:300). Filters were incubated in the dark for 20 min and mounted on glass slides with a drop of 50% phosphate buffer (6.7 mM, pH 7.8) and 50% glycerol, containing 0.5% ascorbic acid. Slides were stored at a temperature of −20 ◦C until analysis. Viral counts were obtained by epifluorescence microscopy (Olympus BX 51, 1250× magnification, equipped with a blue excitation filter, Tokyo, Japan) and were expressed as the number of virus-like particles (VLP).

Bacterial cell production was estimated by measuring the incorporation of 3H-thymidine into bacterial DNA [41]. Methyl-3H-thymidine was added to 10 mL samples at a final concentration of 10 nmol (specific activity: 86 Ci mmol−1). Triplicate samples, together with a formaldehyde-killed adsorption control (final concentration: 0.5%), were incubated for one hour. The incubations were stopped with formaldehyde (final concentration: 0.5%). The thymidine samples were extracted with ice-cold trichloroacetic acid (TCA), according to Fuhrman and Azam [41]. Finally, the TCA-insoluble fraction was collected by filtering the samples through 0.2 μm pore size polycarbonate filters.

The biomass of studied picoplankton groups was estimated using the following cell-to-carbon conversion factors: 20 fgC cell-1 for heterotrophic bacteria [42,43], 36 fgC cell−<sup>1</sup> for *Prochlorococcus* [44], 255 fgC cell−<sup>1</sup> for *Synechococcus* [44], and 2590 fgC cell−<sup>1</sup> for picoeukatyotes [44] and for AAPs [45].

An empirical model was used to examine the regulation of bacteria by predation [46]. Information about coupling between the abundance of bacteria and heterotrophic nanoflagellates (HNF) would be analyzed using a log-log graph. In particular, this graph included Maximum Attainable Abundance (MAA) and Mean Realized Abundance (MRA) lines. The MAA line described the HNF abundance reached at a given bacterial abundance (max log HNF = −2.47 + 1.07 log bacterial abundance). Data close to the MAA line, thus, suggested a strong coupling between the bacteria and HNF abundance, likely interpreted as strong predation on the bacteria [46]. Data positioned below the MRA line instead, suggested that bacterial abundance was not controlled by HNF grazing.

#### *2.3. Statistical Analysis*

Pearson's correlation analysis was carried out to examine the relationship between picoplankton community members and environmental variables. The response of the picoplankton community to environmental conditions was analyzed using multivariate statistical analyses. Principal component analysis (PCA) was performed using CANOCO software (http://www.canoco5.com/), v5 [47]. PCA is a multivariate statistical method mainly used for data reduction. Analysis attempts were made to identify a few components that explain the major variation within data. In our study, PCA was used to better understand the relationships in multivariate data set between environmental and biological parameters in the SAP. The data were centered and standardized to remove the large differences in abundances between groups. Model results were reproduced in ordination biplots, summarizing the main trends in the data.

#### **3. Results**

#### *3.1. Environmental Parameters*

The main physical characteristics observed throughout the water column in December 2015 (Figure 2) and April 2016 (Figure 3) were depicted from the distribution of θ, S, and DO concentration along the transect crossing the SAP (roughly from Bari, Italy toward Dubrovnik, Croatia), and from the vertical profiles of the stations 8 and 9, where picoplankton was sampled. Hereafter, we have referred to the upper layer considering the water column between the surface and 100 m depth, to the intermediate layer considering depths between 100 m and 800 m, and to the deep layer considering depths between 800 m and the bottom (~1250 m depth).

In December 2015, the upper layer along the Ba-Du section was quite heterogeneous, which might be due to the contrasting water masses transported into the SAP by the EAC along its eastern side [18], and by the WAC along its western side (Figure 2). In particular, relatively warm and saline waters moving along the eastern margin of the SAP protruded offshore, reaching the central zone of the pit (Figure 2a,b) due to local cyclonic circulation. There, complex features, such as mesoscale eddies, determine large thermohaline differences among close stations, especially between stations 8 and 9. DO distribution (Figure 2c) in the upper layer showed a marked horizontal gradient, with values diminishing from west to east. The upper intermediate layer, between 100 and 400 m, although more homogeneous, was characterized by the presence of water with properties (θ > 14.30 ◦C and S up to 38.95) typical of the Ionian surface water and the LIWs/CIWs (Levantine/Cretan Intermediate Waters). In the central zone of the SAP, θ gradually decreased with increasing depth, while S had a structure with alternating fresher and saltier layers. Moreover, between 200 and 300 m depth, a branch of fresher water with local S minimum ~38.70 extended from the western flank towards the center of the pit (Figure 2b). In the lower intermediate layer, between 400 and 800 m, instead, S increased up to 38.84, while θ slightly decreased down to 13.60 ◦C. Overall, layers characterized by high S values were also characterized by reduced DO values (Figure 2c). Notably, two relative oxygen minima were found around 100 m and 500 m depth, in correspondence of relative S maxima (Figure 2e,f). However, there were some differences between stations 8 and 9, concerning mostly S and DO. At station 8 between 100 and 300 m, S values were lower than at station 9 in the same layer. The origin of that fresher layer might be attributable to a detachment of the western-coast vein of freshwater originating from the northern/middle Adriatic, as hinted from the S transect. At station 8, there was a DO minimum at 700 m depth, which was not so evident at the nearby station 9. This seems to be connected with the small scale recirculation of the "old" intermediate layers, associated with the LIW/CIW. The deep layer (>800 m) of the SAP was occupied by relatively cold, less saline, and dense waters likely formed during previous winters in the northern Adriatic. These waters would have been trapped in the deepest part of the SAP, which lies below the depth of the Otranto sill [32,34]. Dense waters coming from the northern Adriatic usually sink to the deep part of the SAP following bathymetric constraints over the slope [48] and they reach their equilibrium depths according to their densities. At the time of the cruise

(December 2015), θ and S in the deep layer had values of 13.10 ◦C and 38.71, respectively (Figure 2d,e). As far as the oxygen content is concerned, higher values were found between 800 m and 1000 m depth (Figure 2f), and lower values were found below 1000 m depth. This was a sign of water stagnation and lack of ventilation of the deepest part of the SAP during winters before the oceanographic survey, while, likely, relatively dense waters founded their equilibrium depths between 800 and 1000 m.

**Figure 3.** Vertical distribution of (**a**) potential temperature, (**b**) salinity, (**c**) dissolved oxygen, along the Bari-Dubrovnik (Ba-Du) section during the ESAW2 cruise in April 2016. Grey dots indicate locations of picoplankton sampling. Panels (**d**–**f**) show the vertical profiles of the same properties at stations 8 and 9, except the oxygen data, which were taken from stations 7 and 10 (data from station 8 and 9 are missing, due to technical problems).

In April 2016 (Figure 3), the largest differences with respect to the pre-winter conditions in December 2015 (Figure 2) were found mainly in the upper layer temperature, due to the season signal. The temperature differences between the western and eastern flanks, as well as between the surface and bottom layers over the western shelf, diminished with respect to December 2015 (Figure 3a). On the western slope, at depths between 300 and 500 m (stations 3 and 4), a branch of relatively cold, fresh, and ventilated water (θ ~14.10 ◦C, S ~38.74, DO ~4.9–5.1 mL/L) pointed to a possible intrusion of northern and/or middle Adriatic waters. However, data collected in the deepest layers revealed that θ and DO did not change significantly with respect to the pre-winter period, while S values slightly increased from 38.72 to 38.74 (Figure 3d–f). At the two nearby stations, 8 and 9, the thermohaline properties were almost uniform. From the DO profiles at the stations 7 and 10 (close to 8 and 9, where DO values were not recorded due to technical problems), we might guess that the vertical convection and mixing, resulting in a homogenous vertical distribution of DO, reached 350 m at station 9, and 250 m at station 8. In any case, the vertical distribution of the physical parameters and DO (Figure 3) suggested that a weak vertical convection occurred during winter 2015–2016, and it probably did not exceed 400 m in depth, as confirmed by time-series data recorded at E2-M3A (http://nettuno.ogs.trieste.it/e2-m3a/), also located in the SAP (Figure 1). The highest S values (~38.94) measured in April 2016, associated with the LIW influence, were slightly lower than those observed in December 2015 (~38.95). From the physical parameters, we inferred that the doming structure, typical of the cyclonic circulation in the SAP, was much more enhanced in April than in December. Moreover, it was centered near station 7. This means that the sub-basin-scale cyclonic gyre was probably stronger in April than in December, favoring lateral exchanges along the perimeter of the SAP. The lateral exchange between both coastal flanks and the middle of the transect seemed less active with respect to December 2015. Biogeochemical properties were obtained through chemical analyses performed on water samples collected from discrete depths along the SAP transect. Here we reported the data from stations 8 and 9, relevant for the biological sampling (Figure 1). The vertical distribution of nutrients was similar at the two stations 8 and 9 (Figure 4), for both cruises. Nitrate, phosphates, and silicates were depleted in the upper layer and increased with increasing depth. Highest values of nitrates and phosphates were found between 400 and 600 m, whereas silicates increased almost uniformly until 800 m, and then more rapidly between 1000 m and 1100 m, reaching the highest values in the bottom waters, where turbidity (proportional to the suspended particle concentration) was also high (Figure S1a,c). Ammonia concentrations were higher in December 2015 than in April 2016. Nitrites were low, but local maxima were observed at depths between 50 and 100 m in both December and April, corresponding to fluorescence maxima (Figure S1b,d). The DON and DOP concentrations throughout the water column were higher in December 2015 than in April 2016 (Figure 5). On the contrary, DOC concentrations were higher in April 2016 than in December 2015 (Figure 5a,b). High concentrations of DOC, DON, and DOP were found in the upper layer, between the surface and 100 m. In the upper layer, at 100 m depth, turbidity was relatively high in both periods, with apparently higher values in April 2016 (Figure S1a,c). Turbidity slightly increased also between 200 m and 600 m, where a prominent LIW signal was found.

**Figure 4.** Vertical profiles of nutrients. Station 8: panels (**a**–**e**) and Station 9: panels (**f**–**j**). Periods: December 2015 ESAW1 (black line) and April 2016 ESAW2 (grey line). Symbols indicate sample locations.

**Figure 5.** Vertical profiles of dissolved organic nitrogen (DON), dissolved organic phosphorus (DOP), and dissolved organic carbon (DOC). Station 8: panels (**a**–**c**). Station 9: panels (**d**–**f**). Periods: December 2015 ESAW1 (black line) and April 2016 ESAW2 (grey line). Symbols indicate sample locations.

#### *3.2. Picoplankton Community*

In both surveys, among the measured microbial community (heterotrophic bacteria, *Synechococcus*, *Prochlorococcus*, picoeukaryotes), *Synechococcus* abundances (Table S1) were the highest in the autotrophic fraction and reached values up to 93.93×10<sup>3</sup> cells ML<sup>−</sup>1. *Prochlorococcus* and picoeukaryotes abundances were up to 47.09 <sup>×</sup> <sup>10</sup><sup>3</sup> cells mL−<sup>1</sup> and 2.96 <sup>×</sup> <sup>10</sup><sup>3</sup> cells mL<sup>−</sup>1, respectively. Bacterial abundance ranged from 0.03 <sup>×</sup> <sup>10</sup><sup>6</sup> to 0.24 <sup>×</sup> <sup>10</sup><sup>6</sup> cell mL−<sup>1</sup> with the highest values found above 100 m depth.

In December, the autotrophic biomass (an average of 13.6 μgCL−1) was almost six times higher than heterotrophic (an average of 2.29 μgCL−1), with the domination of *Synechococcus*. Vertical distribution revealed the prevalence of autotrophic biomass over heterotrophic in the epipelagic

layer but also deep waters (Figure 6a,b). In April, the heterotrophic biomass was similar to that in December (on average 2.90 μgCL−1) and slightly higher than the autotrophic biomass (on average 2.47 μgCL<sup>−</sup>1). The autotrophic biomass, mainly composed by picoeukaryotes, was slightly higher than the heterotrophic biomass in the epipelagic layer, while in the deep waters, the heterotrophic biomass dominated among the observed communities (Figure 6c,d).

Bacterial production ranged from 0.01 <sup>×</sup> 104 to 0.09 <sup>×</sup> 10<sup>4</sup> cells h−<sup>1</sup> mL−1. The highest values were found at the surface, but surprisingly high values were recorded at 800 m and 1200 m depth, during April 2016 (Figure 7a). The ratio of nucleic acid content showed the domination of LNA cells in all samples within the bacterial community throughout the euphotic zone (up to 200 m) and the prevalence of HNA cells in the deeper layers (Figure 7b). The abundance of viruses ranged from 0.39 <sup>×</sup> <sup>10</sup><sup>6</sup> VLP mL−<sup>1</sup> to 5.37 <sup>×</sup> 106 VLP mL−<sup>1</sup> (Figure 7c). Overall, their vertical distribution revealed a large abundance in the euphotic zone and low abundance below 300 m depth, at both sampling sites. However, at both stations, values of VLP slightly increased in the layers below 800 m in December, while at station 8, also increased in the deep layer in April. The average virus to bacteria ratio (VBR) was higher in December (19 ± 10 at station 8 and 29 ± 12 at station 9) compared to April (8 ± 2 at station 8 and 8 ± 4 at station 9).

AAPs were observed in all samples during April, with abundance ranging from 0.04 <sup>×</sup> 104 to 0.61 <sup>×</sup> <sup>10</sup><sup>4</sup> cells mL−1. In general, their vertical distribution showed similar patterns at both stations 8 and 9. Below 200 m, in correspondence of the nutricline and the minimum fluorescence layer (Figure S1b,d) the abundances of AAPs started decreasing steeply (Figure 7d). The proportion of AAPs abundances in total prokaryotes ranged from 0.65% to 2.48%, while the proportion of the AAPs biomass ranged from 0.37% to 4.09%, respectively.

**Figure 6.** Vertical profiles of biomass for heterotrophic bacteria (BACT), *Prochlorococcus* (PRO), *Synechococcus* (SYN), picoeukaryotes (PE), aerobic anoxygenic phototrophs (AAP) at stations 8 (**a**) and 9 (**b**) in December 2015 (ESAW1) and at stations 8 (**c**) and 9 (**d**) in April 2016 (ESAW2).

**Figure 7.** Vertical profile of bacterial production (**a**), HNA % (High Nucleic Acid content) (**b**), viruses (VLP) (**c**), and AAPs (**d**) at the station 8 and 9 in December 2015 (ESAW1) and in April 2016 (ESAW2).

#### *3.3. Distribution of Picoplankton Community Members Concerning Environmental Variables*

Relationship between environmental parameters and picoplankton was tested by Pearson's correlation analysis (Table S2). DOC, HNF, and temperature seemed to be the most important environmental factors for bacterial abundance and production, exhibiting significant positive correlations. Ammonium ion displayed a significant positive correlation with *Prochlorococcus* and *Synechococcus.* During the investigation period, VLP significantly correlated with autotrophic members of the picoplankton community.

Positive correlations of AAPs were found for temperature (r = 0.83; *p* < 0.05), for fluorescence (r = 0.76; *p* < 0.05), and for nitrites, nitrates, phosphates, and total dissolved phosphates (r = 0.85; r = −0.94; r = −0.90; r = −0.69; *n* = 14; *p* < 0.05). Furthermore, AAPs were also positively correlated with picoeukaryotes and HNF (r = 0.74; r = 0.91; *n* = 14; *p* < 0.05).

PCA ordination of picoplankton groups concerning environmental variables was used to analyze main factors affecting the abundances of the picoplankton groups (Figure 8a). PCA clustered the samples corresponding to depth (E-epipelagic layer; D-deep layer) and period (Figure 8b). The first principal component explained 44.83% of the variance and correlated positively with HNF and temperature, negatively with a concentration of nitrates and phosphates. The second principal component explained 17.75% of the variance and correlated positively with picoautotrophs and abundance of VLP (Table S3). Picoplankton community members were distributed in two groups, based on the correlation between them. On the one hand, bacteria, bacterial production, and LNA bacteria were closely related, while on the other hand, this was also the case for *Synechococcus* and *Prochlorococcus*. Bacterial cluster showed the strongest correlation with DOC, HNF, and temperature, as indicated by the perpendicular projection of bacterial arrow-tips on the line overlapping the DOC, HNF, and temperature arrow. Cyanobacteria were more related to ammonium, while HNA bacteria showed a positive relationship with nitrates and picoeukaryotes with VLP.

**Figure 8.** PCA (principal component analysis) of (**a**) biological (*Synechococcus*-SYN, *Prochlorococcus*-PRO, picoeukaryotes-PE, total bacteria—HB, HNA HB, LNA HB, VLP, HNF, bacterial production BP) and environmental parameters (Salinity—S, Temperature—T, nutrients). The arrow direction points out to the steepest increase of the variable. The angle between arrows indicates correlations between variables (an angle < 90◦ between two arrows of interest implies positive correlation), whereas the length of an arrow depicts the strength of association between a variable and the ordination axes shown in the biplot. (**b**) PCA clustering of the samples corresponding to depth (E-epipelagic layer; D-deep layer) and period. LNA: Low Nucleic Acid; VLP: virus-like particle; HNF: heterotrophic nanoflagellates.

The relationship between HNF as the main predator of bacteria and the bacteria was described by the Gasol model [46], which we used for our dataset. Most samples from the epipelagic layer were placed above the MRA line, whereas most samples from deep waters were below MRA. This pattern suggested a stronger coupling between bacteria and HNF in the epipelagic layer. In the deep waters layer, HNF predation pressure on bacteria was lower, suggesting larger importance for bottom-up control of bacteria and top-down control of HNF (Figure 9).

**Figure 9.** Relationship between bacterial and HNF abundance at study stations, plotted in a theoretical model (Gasol, 1994) (MAA-maximum attainable abundance: max log HNF = −2.47 + 1.07 log bacterial abundance; MRA-mean realized abundance: mean log HNF = −1.67 + 0.79 log bacterial abundance) in epipelagic (epi-) and deep (deep-) water layers.

#### **4. Discussion**

One of the main findings in this survey was the unusually high biomass of picoautotrophs found in the epipelagic and deep layers in December 2015. *Synechococcus* dominated the picoplankton community with a maximum value of 23.95 μgCL<sup>−</sup>1. Observed high value of *Synechococcus* biomass agrees with those measured in the Levantine basin during the late summer-autumn period [49]. It is well known that picophytoplankton tends to dominate autotrophic biomass and primary production in oligotrophic waters, like those of the Mediterranean Sea [50–52]. Furthermore, picoeukaryotes contributed significantly to the autotrophic biomass due to their larger size and carbon content, compared to cyanobacteria.

Physical and biogeochemical measurements performed in December 2015 revealed the presence of LIW in the intermediate layer and NAdDW at depths > 800 m. Our results also revealed the presence, in the deep layers of the SAP, of *Synechococcus*, *Prochlorococcus,* and picoeukaryotic cells coming from the Eastern Mediterranean and the North Adriatic Sea. Our findings agree well with previous studies [6,53], which associated the high abundances of *Synechococcus* and bacterial cells in the deep layers with water mass movement.

The results of this research showed the remarkably high autotrophic biomass in the deep waters, which is not characteristic of that environment. We suggest that picoautotrophs, aside from their important role in primary production, could be a significant carbon source for higher trophic levels in the form of dead or live prey, or as sinking particles in upper [54–58] and deep ocean waters [59]. In this study, physical and biogeochemical measurements performed in the Southern Adriatic revealed that weak vertical convection occurred during winter 2015–2016, which probably did not exceed 400 m depth. Moreover, a strong heterogeneity of the physical properties distribution, both horizontally and vertically, suggested the large influence of lateral exchanges due to the cyclonic circulation and mesoscale activity, apparently stronger during spring.

In April 2016, the picoplankton community throughout the water column was dominated by heterotrophic biomass (contrary to the December case), which is consistent with previous research for the open Adriatic area [3]. *Synechococcus*, *Prochlorococcus,* and picoeukaryotes abundances fluctuate in a range typical for the euphotic zone in the open Adriatic Sea [3,60,61]. In this survey, *Synechococcus* was the most abundant autotroph. *Synechococcus* hold the advantage over genus *Prochlorococcus* and thrive in P-depleted environments due to the high affinity for inorganic phosphorus and higher phosphate uptake rates as reported recently [62,63]. Furthermore, both genera of cyanobacteria showed a positive relationship with the concentration of ammonium, which is supported by the fact that *Synechococcus* can exploit both the reduced and the oxidized forms of nitrogen. In oligotrophic environments, ammonium represents the major nitrogen compound for primary production [64], and *Prochlorococcus* are effective in consuming the same [65,66].

Our findings clearly showed that bacterial production and biomass values measured in the deep waters were similar to the values for the epipelagic layer. La Ferla et al. [67] also found similar vertical patterns of bacterial density and bacterial production in the Western and the Eastern Mediterranean Sea. Additionally, we determined the dominance of HNA bacterial group in the deep-sea, while the LNA bacterial group prevailed in the epipelagic layer. Recent research from the Adriatic Sea pointed out the increase of HNA bacteria with depth, and together with LNA bacteria, they participated in the bacterial activity [10,68]. It is well known that HNA bacterial cells are larger than LNA cells and that the proportion of HNA bacteria increases with depth [69,70]. Larger cell size is reflected in a higher-level cell-specific activity or higher metabolic versatility in the mesopelagic ocean [5,70,71]. So, we argue that similar values of bacterial production measured in the deep water and epipelagic layer, together with the dominance of HNA bacteria in deep water, showed bacterial activity in the deep waters of SAP. We can conclude that the results of bacterial production as the measure of activity in this study support the idea that deep ocean prokaryotes are active as those living in the epipelagic waters [4,72–75]. Beside the domination of HNA cells in the bacterial community, our data showed increased values of DOC, phosphates, and bacterial production below 800 m during the post-winter period. This suggests

that the bacterial activity increases with the supply of organic carbon, as Nagata et al. [76] previously described for the deep-sea environment.

Bacterial abundance throughout the water column was approximately 10<sup>5</sup> cells mL<sup>−</sup>1, which is per the recent data reported for the SAP [8,10] and with previous studies carried out in the Mediterranean Sea [67,72,77–79]. Heterotrophic bacteria contribute at a high percentage in total picoplankton biomass, acting as decomposers of organic matter and important producers of new biomass. Through grazing by HNF or larger zooplankton, their biomass becomes available to higher trophic levels, in both the epipelagic and deep waters of the Mediterranean Sea [5,46,59]. During our surveys, bacterial abundance and productivity showed positive relationships with DOC concentration and HNF abundance. Our results showed that the increase in bacterial abundance and cell production supported the increase in the number of HNF, especially in the epipelagic layer. It revealed that bacteria constitute a potential food resource for the nanoflagellate community and suggest a strong top-down control of bacteria. These results confirm previous findings showing [2,80] that predators prefer active bacteria and remove bacterial production, and that they can control the abundance of the bacterial community in surface waters [46]. Some studies [21,81] suggested strong bottom-up control of bacteria, which is in accordance with our results, showing a positive relationship between bacterial parameters and DOC concentration through the water column.

Little is known about the role and distribution of AAPs in the Adriatic Sea, especially in its deepest part. Šanti´c et al. [12] investigated AAPs abundances along the eastern Adriatic, in coastal and transitional waters, and found that their proportion in total prokaryotes was 7.3% (±4.3%). Celussi et al. [11] pointed out that the concentration of bacteriochlorophyll *a*, the main pigment of AAPs, could be converted into abundance. Therefore, AAPs might represent up to 10% of total prokaryotes in the open Adriatic Sea. In recent research of AAPs, the abundance, the proportion in total prokaryotes, and biomass decreased along the trophic gradient [45]. Furthermore, these counts were substantially low at the oligotrophic open sea station, especially under 70 m depth. Similar vertical distribution was already noted [13,82–86] since these organisms inhabit the euphotic zone because of their phototrophic nature. The results of this study are consistent with previous reports, which also found lower abundances in the more oligotrophic area when compared to the more productive regions or shelf seas [13,39,82,85–89].

The main environmental parameters that influenced the picoplankton community according to the research conducted by Šanti´c et al. [12] were chlorophyll, followed by nitrates and temperature.

In this research, fluorescence, as a proxy of chlorophyll, nitrates, nitrites, phosphates, and total dissolved phosphates, influenced the abundance of AAPs. Besides, the correlation found between bacteria and AAPs, which has also been reported by Celussi et al. [11], could be explained by the fact that both AAPs and heterotrophic bacteria rely on labile DOC, part of which is released by phytoplankton. Secondly, AAPs and phytoplankton have a similar dependence on light [90], and this relationship may reflect the same dependence on limiting nutrients, such as phosphorus or nitrogen [15]. The majority of AAP surveys have shown that temperature, salinity, chlorophyll concentration, and carbon and nitrogen compounds affect the dynamics and distribution of the AAP population [39,85–87,89]. The high positive correlation between AAPs and HNF found during our study suggest that AAPs represent an important prey for flagellates in the open Adriatic Sea, as documented in the Mediterranean Sea [91,92]. It is a well-known fact that AAPs are larger than the average heterotrophic bacteria [13,14,86,93] and that, consequently, their biomass contribution to total bacterial biomass is also high (up to 11% in this survey). Cell size is an important parameter in feeding ecology, which provides a valuable insight into the trophic significance of AAPs. It was previously observed that nanoflagellates preferentially ingest the larger bacterial cells [94–96].

The number of viruses ranged between 0.39 <sup>×</sup> 10<sup>6</sup> VLP Ml−<sup>1</sup> and 5.37 <sup>×</sup> 10<sup>6</sup> VLP mL<sup>−</sup>1, which is in accordance with other studies for the Adriatic and Mediterranean Sea [97–101]. Their average abundance was 16 ± 12× higher than the abundance of prokaryotes, as well as in previous studies on oligotrophic deep waters [99]. Relatively high VBR determined in the deep waters could be a result of higher lytic activity caused by lysogenic viruses. Water mass sinking to the SAP could cause the induction of lysogenic viruses and thus keep the viral abundance stable. Winter et al. [102] recently pointed out that mixing in the deeper waters led to the induction of lysogenic viruses. Viruses positively correlated with the members of the picoplankton community at the sampling sites, indicating their involvement in shaping the picoplankton abundances in the oligotrophic waters of the SAP. The PCA analysis revealed a close relationship between viruses and the picoautotrophic group. It is known that viruses can be an important mortality factor for picoeukaryotes in oligotrophic oceanic waters [103]. Our results suggest that viruses are an important component of the microbial food web, but this topic should be addressed in further studies.

#### **5. Conclusions**

The present study points out the importance of *Synechococcus*, *Prochlorococcus*, picoeukaryotes, heterotrophic bacteria, AAPs, heterotrophic nanoflagellates, and viruses altogether in the carbon flux through the microbial food web in the SAP for the first time, which depended largely on water mass movements and mixing.

A significant amount of picoautotropic biomass was found in the epipelagic, as well as in the deep layers, during the period under the influence of LIW inflow and NAdDW sinking, respectively. After weak vertical convection during winter, heterotrophic bacteria dominated picoplankton biomass throughout the water column. Furthermore, bacterial activity values in the deep layer were similar to the ones obtained in the epipelagic layer. Besides the physical pump, the changes in the picoplankton community in the SAP were mainly driven by nutrient availability, more specifically bacteria by DOC and cyanobacteria by ammonium ion. Our results also reveal that bacteria represent an important prey for nanoflagellates, especially in the epipelagic layer. Moreover, autotrophic and heterotrophic members of the picoplankton community could be available prey and a source of carbon for the food web below the epipelagic layer. Finally, AAPs contributed up to 4% of bacterial biomass, which can be transferred through higher trophic levels by grazing. Viral lysis may affect the activity of the picoplankton community and serve as an important source of DOC in the deep-sea.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4441/11/8/1655/s1, Figure S1:Vertical profiles of optical backscattering and fluorescence. Station 8: panels (a) and (b). Station 9: panels (c) and (d). Periods: December 2015 ESAW1 (black line) and April 2016 ESAW2 (grey line), Table S1: Summary of the picoplankton community abundances, Table S2. Summary of the correlations between picoplankton community (Synechococcus-SYN, Prochlorococcus-PRO, picoeukaryotes-PE, total bacteria HB, HNA HB, LNA HB, VLP, HNF, bacterial production BP) and environmental variables (Temperature-T, nutrients). Acronym n.s. stands for not significant, Table S3. Summary of the principal component analysis (PCA) of the picoplankton groups and environmental variables.

**Author Contributions:** Conceptualization, M.B.; Data curation, V.K.; Formal analysis, D.Š. and A.V.T.; Funding acquisition, V.K.; Investigation, D.Š.; Methodology, D.Š., M.G., A.V.T., M.O., C.S., S.Š. and B.G.; Project administration, V.K.; Supervision, M.Š.; Writing—original draft, D.Š.; Writing—review and editing, V.K., M.B., M.G., A.V.T. and M.O.

**Funding:** The study was carried out within the framework of the EUROFLEETS2 research infrastructures project under the 7th Framework Program of the European Commission (grant agreement 312762). It was partially funded by the Italian national project RITMARE (Ricerca Italiana per il Mare, grant numbers: SP3-WP3-AZ1; SP4-LI4-WP1; SP5-WP3-AZ3) and by the Croatian Science Foundation as a part of research projects: IP-2014-09-4143 "Marine microbial food web processes in global warming perspective" (MICROGLOB).

**Acknowledgments:** We thank L.U. for nutrients, TDN and TDP analyses, H.M. and S.M. for their help with sampling, and the crew of R/V BIOS DVA.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Numerical Investigations of Tsunami Run-Up and Flow Structure on Coastal Vegetated Beaches**

#### **Hongxing Zhang 1, Mingliang Zhang 1,\*, Tianping Xu <sup>1</sup> and Jun Tang <sup>2</sup>**


Received: 23 October 2018; Accepted: 29 November 2018; Published: 3 December 2018

**Abstract:** Tsunami waves become hazardous when they reach the coast. In South and Southeast Asian countries, coastal forest is widely utilized as a natural approach to mitigate tsunami damage. In this study, a depth-integrated numerical model was established to simulate wave propagation in a coastal region with and without forest cover. This numerical model was based on a finite volume Roe-type scheme, and was developed to solve the governing equations with the option of treating either a wet or dry wave front boundary. The governing equations were modified by adding a drag force term caused by vegetation. First, the model was validated for the case of solitary wave (breaking and non-breaking) run-up and run-down on a sloping beach, and long periodic wave propagation was investigated on a partially vegetated beach. The simulated results agree well with the measured data. Further, tsunami wave propagation on an actual-scale slope covered by coastal forest *Pandanus odoratissimus* (*P. odoratissimus*) and *Casuarina equisetifolia* (*C. equisetifolia*) was simulated to elucidate the influence of vegetation on tsunami mitigation with a different forest open gap. The numerical results revealed that coastal vegetation on sloping beach has significant potential to mitigate the impacts from tsunami waves by acting as a buffer zone. Coastal vegetation with open gaps causes the peak flow velocity at the exit of the gap to increase, and reduces the peak flow velocity behind the forest. Compared to a forest with open gaps in a linear arrangement, specific arrangements of gaps in the forest can increase the energy attenuation from tsunami wave. The results also showed that different cost-effective natural strategies in varying forest parameters including vegetation collocations, densities, and growth stages had significant impacts in reducing the severity of tsunami damage.

**Keywords:** tsunami waves; numerical simulation; wave run-up; flow structure; coastal vegetation

#### **1. Introduction**

Tsunamis are generated by marine earthquakes, underwater volcanic eruptions, or submerged landsides. Recent tsunami disasters (e.g., in the Indian Ocean in 2004, Samoa in 2009, Chile in 2010, Tohoku in 2011, and Indonesia in 2018) have caused vast losses of life and property, destruction of critical infrastructure in low-lying coastal areas, and massive damage to the coastal ecosystem [1–4]. These great threats demonstrate the need to mitigate tsunamis by using artificial obstacles including wave-dissipating blocks, rock breakwaters, and large embankments. Recent studies investigated tsunami waves in both laboratory experiments [5–8] and field investigations [9,10] with the ultimate goal of mitigating their impacts and reducing natural disasters in coastal regions. Considering the advantages of flexibility and repeatability in numerical models, many numerical simulations have been carried out using different models to investigate the run-up heights and velocity processes of tsunamis [11–16]. Typical models that were considered included depth-integrated Boussinesq-type

models and depth-integrated non-linear shallow water models (NSWMs), with results showing that NSWMs are more robust and efficient for tsunami predictions in real-scale basins.

Generally, coastal forests which act as natural obstacles are widely distributed in coastal areas. They also effectively mitigate tsunami damage from economic, environmental, and aesthetically pleasing points of view in concert with hard structures [17,18]. Currently, there is some literature on the capability of coastal vegetation to attenuate wave energy in short-period waves [19,20]. While tsunamis are categorically different from short-period waves, there have been an increasing number of studies that address the role of vegetation in mitigating coastal natural disasters, which can include the strategy of planting coastal vegetation as a bio-shield [21–23]. Recent research has calculated the friction and drag force of vegetation in a tsunami, using mathematical equations. Real-scale simulations of tsunami were performed to investigate the effectiveness of forests as a bio-shield for tsunami protection, with results indicating that a forest with two layers in the vertical direction including *P. odoratissimus* and *C. equisetifolia* was effective for attenuating the wave energy [24]. Based on two dimensional nonlinear long-wave equations, Thuy et al. developed a numerical model to estimate tree-breaking, the drag forces, and turbulence-induced shear force due to the presence of vegetation (*P. odoratissimus* and *C. equisetifolia*) [25]. Considering the effects of the porosity of vegetation, a one-dimensional numerical model using Boussinesq-type equations was developed to evaluate the effects of the forest density distribution on tsunami-force reduction, including the drag and inertia forces caused by vegetation [26]. Numerical investigations were conducted to study the impact of patchy vegetation on tsunami dynamics based on Boussinesq model, and the results demonstrated that patchy vegetation, with appropriate configuration, can be effective in mitigating tsunami hazards [27]. Based on a nonlinear long wave equation model that included the breaking or washout of trees, numerical simulations were carried out to estimate the effects of coastal forest and sea embankments on reducing the washout area and the tsunami mitigation function of the coastal forest [28].

Roads are perpendicular to the coast function as pathways that form gaps in the vegetation of the coastal forest. Gaps in the coastal forest can increase risks and potential damage, as the water flow from the tsunami accelerates as it moves through the gap into the densely populated block [29]. Tanimoto et al. studied it by changing the width of open gap and found that a specific gap with 15 m width causes the highest flow velocity in their simulated conditions [30]. Thuy et al. numerically simulated the effect of an open gap in coastal forests on tsunami run-up, and found that maximum flow velocity greatly increased at the open gap exit, meaning that an open gap (like a road) in a coastal forest had a negative effect on tsunami run-up behind the forest [31]. Nandasena et al. investigated hydrodynamic parameters (flow velocity and depth) of tsunami waves and bending moment of vegetation on Misawa, a site covered by pine forest with two gaps [32]. However, no studies have discussed the effects of different open gaps in existing forests via the analysis of maximum tsunami run-up and variations of velocity through the gaps using a model with shock capture capability.

In this paper, we numerically investigate the hydrodynamic processes on a vegetated sloping beach and quantify the effects of a gap in forest vegetation on the run-up by solving the depth-averaged 2D model. The proposed model was tested for breaking and non-breaking solitary waves propagating on a bare sloping beach, and long periodic wave propagation on a partially-vegetated sloping beach to examine the accuracy of the numerical model. Then, the model was used to simulate tsunami waves propagating on actual-scale forest-covered beach, and to study the mitigation effects of gap arrangements and different vegetation parameters on tsunami waves.

#### **2. Numerical Method**

#### *2.1. Governing Equations*

The depth-averaged 2D shallow water equations are formed by integrating the Navier–Stokes equations which include continuity and momentum equations for depth-averaged free surface flows. These are expressed as Equations (1)–(3):

$$\frac{\partial \mathbf{h}}{\partial t} + \frac{\partial uh}{\partial \mathbf{x}} + \frac{\partial vh}{\partial y} = 0 \tag{1}$$

$$\frac{\partial \text{h}\mathbf{u}}{\partial t} + \frac{\partial (\text{h}\mathbf{u}\mathbf{u})}{\partial \mathbf{x}} + \frac{\partial (\text{h}\mathbf{v}\mathbf{v})}{\partial y} - \frac{\partial}{\partial \mathbf{x}}(v\_l h \frac{\partial \mathbf{u}}{\partial \mathbf{x}}) - \frac{\partial}{\partial y}(v\_l h \frac{\partial \mathbf{u}}{\partial y}) = -g h \frac{\partial \eta}{\partial \mathbf{x}} - \tau\_{\mathbf{b}\mathbf{x}} + f\_c l v - f\_{\mathbf{x}} \tag{2}$$

$$\frac{\partial lv}{\partial t} + \frac{\partial (hvv)}{\partial x} + \frac{\partial (hvv)}{\partial y} - \frac{\partial}{\partial x}(v\_lh\frac{\partial v}{\partial x}) - \frac{\partial}{\partial y}(v\_lh\frac{\partial v}{\partial y}) = -gh\frac{\partial \eta}{\partial y} - \tau\_{by} - f\_c hu - f\_y \tag{3}$$

In the above equations, *t* is time, *h* indicates the local water depth, *u* and *v* stand for the depth-averaged flow velocities in the *x* and *y* directions, respectively, *η* stands for the water surface elevation from a reference datum, *vt* is defined as the eddy viscosity coefficient calculated by *ν<sup>t</sup>* = *αu*∗*h* where α is empirical constant and ranges from 0.3 to 1.0, *u*<sup>∗</sup> means the bed shear velocity, *τbx* and *τby* denote friction in the *<sup>x</sup>* and *<sup>y</sup>* directions, respectively, *<sup>τ</sup>bx* <sup>=</sup> <sup>g</sup>*n*2*<sup>u</sup>* √ *u*<sup>2</sup> + *v*<sup>2</sup> *h* 1 3 , *<sup>τ</sup>by* <sup>=</sup> <sup>g</sup>*n*2*<sup>v</sup>* √ *u*<sup>2</sup> + *v*<sup>2</sup> *h* 1 3 , where *n* stands for Manning's roughness coefficient, *fc* indicates the Coriolis parameter, and *fx* and *f* y stand for the drag force caused by vegetation.

The conservative and vector form of 2D shallow water equations with vegetation effects are expressed as follows:

$$\frac{\partial \mathbf{U}}{\partial t} + \frac{\partial \mathbf{F}}{\partial \mathbf{x}} + \frac{\partial \mathbf{G}}{\partial y} = \frac{\partial \mathbf{F}\_d}{\partial \mathbf{x}} + \frac{\partial \mathbf{G}\_d}{\partial y} + \mathbf{S} \tag{4}$$

The definitions of **U**, **F**, **G**, **F***<sup>d</sup>* and **G***<sup>d</sup>* are shown as follows:

$$\mathbf{U} = \begin{bmatrix} h \\ hu \\ hv \end{bmatrix}, \mathbf{F} = \begin{bmatrix} hu \\ huu \\ huv \end{bmatrix}, \mathbf{G} = \begin{bmatrix} hv \\ huv \\ hvv \end{bmatrix}, \mathbf{F}\_d = \begin{bmatrix} 0 \\ \nu\_l \frac{\partial u}{\partial x} \\ \nu\_l \frac{\partial u}{\partial x} \end{bmatrix}, \mathbf{G}\_d = \begin{bmatrix} 0 \\ \nu\_l \frac{\partial u}{\partial y} \\ \nu\_l \frac{\partial u}{\partial y} \end{bmatrix}, \mathbf{S} = \begin{bmatrix} 0 \\ -gh\frac{\partial v}{\partial x} - \eta\_{ll} + f\_l hv - f\_x \\ -gh\frac{\partial v}{\partial y} - \eta\_{yy} - f\_t hu - f\_y \end{bmatrix} \tag{5}$$

For notational convenience, Equation (4) is often rewritten as:

$$\frac{\partial \mathbf{U}}{\partial t} + \nabla \cdot \mathbf{E}\_w^{adv} = \nabla \cdot \mathbf{E}\_w^{diff} + \mathbf{S} \tag{6}$$

where **E***adv <sup>w</sup>* <sup>=</sup> **<sup>F</sup>***<sup>i</sup>* <sup>+</sup> **<sup>G</sup>***<sup>j</sup>* and **<sup>E</sup>***di f <sup>w</sup>* = **F***di* + **G***<sup>d</sup> j*.

#### *2.2. Vegetation Drag Force*

The resistance effects of vegetation on flow are added into the momentum equations as an internal source of resistant force per unit fluid mass. The drag force exerted by vegetation per unit volume is expressed as [33]:

$$f\_x = \frac{1}{2} \text{NC}\_D(h) b \text{vmin}(h\_\upsilon, h) u \sqrt{u^2 + v^2}, \\ f\_y = \frac{1}{2} \text{NC}\_D(\mathbf{h}) b\_\upsilon \text{min}(h\_\upsilon, h) v \sqrt{u^2 + v^2} \tag{7}$$

Here, *N* stands for the vegetation density defined as the number of plants per square meter, *bv* indicates the diameter, and *hv* indicates the height of the vegetation. *CD (h)* indicates the

depth-averaged equivalent drag coefficient which considers the vertical stand structures of a tree, defined by Tanaka et al. [1] as follows:

$$\mathbb{C}\_{D}(h) = \mathbb{C}\_{D-r\varepsilon f} \frac{1}{h} \int\_{0}^{h} a(z\_{\mathcal{S}}) \beta(z\_{\mathcal{S}}) dz\_{\mathcal{S}} \tag{8}$$

$$\alpha(z\_{\mathcal{S}}) = \frac{b\left(z\_{\mathcal{S}}\right)}{b\_{rcf}} \tag{9}$$

$$\beta(z\_{\S}) = \frac{\mathbb{C}\_D(z\_{\S})}{\mathbb{C}\_{D \cdot ref}} \tag{10}$$

where *CD-ref* is the reference drag coefficient of the trunk at *zg* (equal to 1.2 m in principle), *α*(zg) is considered as an additional coefficient to express the effects of cumulative width on drag force at each height *zg*, *b*(*zg*) means the projected width, and *bref* indicates the reference projected width. *β*(*zg*) is expressed as an additional coefficient representing the effect of leaves or aerial roots on drag force, and *CD*(*zg*) is drag coefficient of a tree at the height *zg* above the ground surface.

#### *2.3. Finite Volume Method*

The discretization of the governing equations is based on the finite volume method using an unstructured triangular mesh. The conserved variables are defined at the cell centers and represent the average value of each cell. The integral form of Equation (10) over the *i*th control volume can be expressed as:

$$\int\_{V\_i} \frac{\partial \mathbf{U}}{\partial t} d\_V + \int\_{V\_i} \nabla \cdot \mathbf{E}\_w^{adv} d\_V = \int\_{V\_i} \nabla \cdot \mathbf{E}\_w^{diff} d\_V + \int\_{V\_i} \mathbf{S} d\_V \tag{11}$$

In this expression, subscript *i* and *j* denote the element and the element side, respectively, while *Vi* stands for the domain of the *i*th. Basing Green's theorem. Equation (11) can be rewritten as:

$$\frac{\Delta \mathbf{U}\_i}{\Delta t} A\_i = -\oint\_{L\_i} \mathbf{E}\_w^{adv} \cdot \mathbf{n} dl \, + \oint\_{L\_i} \mathbf{E}\_w^{diff} \cdot \mathbf{n} dl \, + \int\_{V\_i} \mathbf{S} d\_V \tag{12}$$

where U*<sup>i</sup>* stands for the average value of the conserved variables over the *i*th cell, and U*<sup>i</sup>* is stored at the center of the *i*th cell, with **U***<sup>i</sup>* = <sup>1</sup> *Ai Vi* **U***dV*. *Li* means the boundary of the *Vi*. *Ai* denotes the area of the *i*th cell, **n** is the outward surface normal vector of *Li*, with **n** = (*nx*, *ny*)=(cos *φ*, sin *φ*), and *φ* is the angle included between the *x* direction and the outward normal vector.

In the 2D triangular grid system, the line integral term in the Equation (12) can be further approximated and assessed as follows:

$$
\Delta \mathbf{U}\_i = -\frac{\Delta t}{A\_i} \sum\_{j=1}^m \left( \mathbf{E}\_w^{adv} \cdot \mathbf{n}\_{i\bar{j}} \right) l\_{i\bar{j}} + \ \ \ \frac{\Delta t}{A\_i} \sum\_{j=1}^m \left( \mathbf{E}\_w^{diff} \cdot \mathbf{n}\_{i\bar{j}} \right) l\_{i\bar{j}} + \ \ \ \frac{\Delta t}{A\_i} \int\_{\tilde{V}\_i} \mathbf{S} d\_V \tag{13}
$$

where *m* denotes the number of total edges of the triangular cell (three in this model), the subscript *j* means the index of the edge of a triangular mesh, **n***ij* stands for the outward normal flux vector, and *lij* means the length of the arc.

#### *2.4. Evaluation of Numerical Fluxes*

Variables are usually approximated as constant states within each control volume. Riemann problems at the interface of the cell can be solved using various Riemann approximations for evaluating the interface fluxes. The interface fluxes of Roe's solver are expressed as follows:

$$\mathbf{E}\_w^{adv} \cdot \mathbf{n} = \frac{1}{2} [ (\mathbf{F}, \mathbf{G})\_R \cdot \mathbf{n} + (\mathbf{F}, \mathbf{G})\_L \cdot \mathbf{n} - |I| (\mathbf{U}\_R - \mathbf{U}\_L) ] \tag{14}$$

where **U***<sup>R</sup>* and **U***<sup>L</sup>* are reconstructed Riemann state variables on the right and left sides of the cell interface, respectively. The flux Jacobian matrix *A* can be assessed as:

$$|I| = \frac{\partial \left( \mathbf{E}\_x^{\text{adv}} \cdot \mathbf{n} \right)}{\partial \mathbf{U}} = \frac{\partial \mathbf{E}}{\partial \mathbf{U}} n\_x + \frac{\partial \mathbf{G}}{\partial \mathbf{U}} n\_y = \begin{bmatrix} 0 & n\_x & n\_y \\ \ \left( \mathbf{c}^2 - \mathbf{u}^2 \right) n\_x - \mathbf{u} \mathbf{v} n\_y & 2 \boldsymbol{u} n\_x + \mathbf{v} n\_y & \boldsymbol{u} n\_y \\ -\mathbf{u} \boldsymbol{v} n\_x + \left( \mathbf{c}^2 - \mathbf{v}^2 \right) n\_y & \boldsymbol{v} n\_x & \boldsymbol{u} n\_x + 2 \boldsymbol{v} n\_y \end{bmatrix} \tag{15}$$

where *nx* and *ny* denote the components of the outward surface normal vector in the *x-* and *y-* directions, respectively, and *c* indicates the wave velocity, with

$$c = \sqrt{\mathfrak{g}^{\text{fr}}}$$

where |*J*| = *R*|Λ|*L*, *R* and *L* stand for the right and left eigenvector matrices, and |Λ| denotes the diagonal matrix of the absolute values of the eigenvalues of *A*.

Where |Λ| is defined as:

$$|\Lambda| = \begin{pmatrix} \Lambda^1 & 0 & 0 \\ 0 & \Lambda^2 & 0 \\ 0 & 0 & \Lambda^3 \end{pmatrix}$$

where

$$
\lambda^1 = \iota m\_x + \upsilon \iota\_{y\_{\mathcal{I}}} \; \lambda^2 = \iota m\_x + \upsilon \iota\_y - \mathfrak{c}, \; \lambda^3 = \iota m\_x + \upsilon \iota\_y + \mathcal{c} \tag{16}
$$

The right and left eigenvector matrices are expressed as follows:

$$R = \begin{pmatrix} 0 & 1 & 1 \\ n\_y & u - c n\_x & u + c n\_x \\ -n\_x & v - c n\_y & v - c n\_y \end{pmatrix}, \\ L = \begin{pmatrix} -(\imath n\_y - \imath n\_x) & n\_y & -n\_x \\ \imath n\_x + \imath n\_y & +\frac{1}{2} & \frac{-n\_x}{2\varepsilon} & \frac{-n\_y}{2\varepsilon} \\ -\frac{\imath n\_x + \imath n\_y}{2\varepsilon} & +\frac{1}{2} & \frac{n\_x}{2\varepsilon} & \frac{n\_y}{2\varepsilon} \end{pmatrix} \tag{17}$$

where the Riemann state variables *u*, *v*, and *c* on the cell interface are necessary to deal with the fluxes as calculated by Roe's average:

$$u = \frac{\sqrt{\hbar\_+} u\_+ + \sqrt{\hbar\_-} u\_-}{\sqrt{\hbar\_+} + \sqrt{\hbar\_-}}, \ v = \frac{\sqrt{\hbar\_+} v\_+ + \sqrt{\hbar\_-} v\_-}{\sqrt{\hbar\_+} + \sqrt{\hbar\_-}}, \ c = \sqrt{\frac{\mathcal{g}(\hbar\_+ + \hbar\_-)}{2}} \tag{18}$$

The subscripts + and − indicate the right and left sides of the cell edge, respectively. When the drying-wetting interfaces exist in the computational domain, Roe's average can be calculated by: *u* = *<sup>u</sup>*<sup>+</sup> <sup>+</sup> *<sup>u</sup>*<sup>−</sup> <sup>2</sup> , *<sup>v</sup>* <sup>=</sup> *<sup>v</sup>*<sup>+</sup> <sup>+</sup> *<sup>v</sup>*<sup>−</sup> <sup>2</sup> .

#### *2.5. Treatment of Wetting and Drying Fronts*

A technique to treat both wet and dry boundaries was introduced to achieve zero mass error [34,35]. A criterion, ε, was adopted to define and classify the following four types of edges:


According to these four types of edges, cells were correspondingly divided into three types:


**Figure 1.** Schematic diagram of wet-dry edges.

#### **3. Numerical Simulation and Experimental Validation**

#### *3.1. Solitary Wave Run-up on a Bare Sloping Beach*

The depth-integrated shallow water model was run for experiment cases of breaking and non-breaking solitary waves on a bare sloping beach conducted by Synolakis to validate the accuracy of the numerical scheme for modeling wave run-up and run-down [36]. The topography consists of a 1:19.85 sloping beach adjacent to a constant depth region, as illustrated in Figure 2. A solitary wave propagates from left to right in a wave flume; and according to first-order solitary wave theory, the water surface elevation and velocity in initial time are defined as follows:

$$\eta(\mathbf{x}) = H\_{\text{Wsec}} h^2 \left( \sqrt{\frac{3HW}{4h\_0^3}} (\mathbf{x} - \mathbf{X}\_0) \right) \tag{19}$$

$$
\mu(\mathbf{x}) = \eta(\mathbf{x}) \sqrt{\frac{\mathcal{S}}{h\_0}} \tag{20}
$$

where *Hw* means wave height, *h*<sup>0</sup> denotes initial water depth with a value of 1 m, *X*<sup>0</sup> stands for the position of initial wave crest and is located a half wave length from the toe of the sloping beach in the computing domain, and *u* is wave velocity.

**Figure 2.** Sketch of the solitary wave propagation on a sloping beach.

In this case study, the numerical model used uniform triangular cells with a grid spacing of 0.02 m and a time step of 0.001 s. The minimum water depth to define a "dry bed" was set as 0.0001 m, and the Manning's bed roughness coefficient was calibrated as 0.01 for the glassed flume in the present model. The solitary wave was defined as non-breaking in this case, as *H*w/*h*<sup>0</sup> = 0.0185. For convenience, the results are presented in non-dimensional forms: x<sup>∗</sup> = *<sup>x</sup> h*0 , *η*<sup>∗</sup> = *<sup>η</sup> <sup>h</sup>*<sup>0</sup> and t ∗ = *t* g *h*0 . Comparisons of the simulated and experimental free-surface evolutions are presented in Figure 3. As illustrated, the incident wave propagates on the sloping beach at the early stage (t∗ = 25, 30, 35, 40, 45, and 50), and reaches maximum run-up height at about t∗ = 55, at which point backwash occurs. The maximum run-down happens at around t∗ = 70. Simulated free surface profiles show good agreement with experimental data. Figure 4 shows the comparison of simulated and experimental water surface processes from a breaking solitary wave with *Hw/h*<sup>0</sup> = 0.3. In Figure 4, wave breaking is not well reproduced by the depth-averaged shallow model and the computed wave fronts are steeper and slightly earlier than the experimental results at t∗ = 15 and t∗ = 20, as the model does not consider a wave dispersion term [37,38]. However, the wave breaking is only in a small portion of the domain. At the next step (t∗ = 25), wave breaking is reasonably simulated by the model as a collapse of the wave near the shoreline approximately which can be explained that the wave breaking may be weakened and become smaller portion of the domain. The breaking solitary wave reaches a maximum height around t ∗ = 45, and approaches the lowest position where a hydraulic jump is formed near the shoreline around t∗ = 55. The simulated results agree with experimental data very well, which indicate that the proposed model can accurately predict the propagation of breaking and non-breaking solitary waves on a bare sloping beach.

**Figure 3.** Run-up and run-down of *Hw/h*<sup>0</sup> = 0.0185 non-breaking solitary wave on a 1:19.85 sloping beach.

*Water* **2018**, *10*, 1776

**Figure 4.** Run-up and run-down of *Hw/h*<sup>0</sup> = 0.3 breaking solitary wave on a 1:19.85 sloping beach.

#### *3.2. Propagation of Long Periodic Waves on a Partially Vegetated Sloping Beach*

A laboratory experiment exploring the propagation of long periodic waves in a wave channel with a partially-vegetated sloping beach was conducted in Saitama University [31,32]. In this experiment, the wave tank was 15 m long and 0.4 m wide for propagation, and the bed profile of the beach consisted of several sections with various slopes (Figure 5a). The vegetation was modelled using wooden cylinders with a diameter of 0.005 m at a still water depth of 0.44 m. The vegetation domain was from *x* = 10.36 to 11.36 m, with different vegetation gap widths (*B*g = 0, 0.07, and 0.4 m), while the vegetation density was set as 2200 stem/m−2, and the drag coefficient *CD* of wooden cylinders was calibrated as 2.5 in the present study. An incident sinusoidal wave with a period of *T* = 20 s and wave height of 0.16 cm was propagated from a 0.52 m long constant bottom segment to a sloping beach. In the current case, water surface elevations were measured using capacitance wave gauges at six positions (G1–G6) along the center of the flume in cases of *B*<sup>g</sup> = 0 and 0.4 m, and the flow velocities were measured using electromagnetic current meters at locations in the cross section passing behind vegetation domain G6 (see Figure 5b) in case of *B*<sup>g</sup> = 0.07 m. The run-up height above still water surface was measured by tracing the moving of water front by eyes with a scale on the slope. The computational domain was represented by uniform triangular cells with a grid spacing of 0.01 m and a time step of 0.002 s for the numerical model. The minimum depth criterion for a dry and wet bed was considered as 0.001 m, and the Manning coefficient was set as 0.012.

**Figure 5.** Experimental setup of wave flume. (**a**) Longitudinal section, (**b**) plan view of vegetation zone and measurement points.

A comparison of wave crests, heights, and troughs in Figure 6a,b demonstrates that for the case of non-vegetated areas, wave heights increased due to shoaling as the long waves approached shallow water, however, the increased amplitudes of wave heights were attenuated by the resistance effects of vegetation at the sloping zone, with the attenuated rate reaching 38.96% at G6 (at the center behind vegetation). The wave height at the front of the vegetated zone was also affected by the reflection

of waves caused by the vegetation. Figure 7 shows the comparison of time series of velocities at the center of the gap exit and the center of the end of the vegetation when *B*<sup>g</sup> = 0.07 m, indicating that the peak flow velocity at the center of the gap exit can reach 0.42 m/s and is 3.07 times more than the peak flow velocity at behind the center of the vegetation zone. Figure 8 shows the distribution of the peak velocities averaged from five wave periods at steady state in a cross-section of Gage 6, indicating that vegetation plays an important role in the distribution of flow velocities. The results demonstrate that the present model is an effective tool to predict long periodic wave propagation on a partially vegetated sloping beach. Vegetation can effectively attenuate the wave propagation; however, an open gap in vegetation zone generates large flow and adverse effect at gap exit.

(**b**)

**Figure 7.** Time series of velocities, (**a**) at the center of the gap exit, (**b**) at the center of the end of the vegetation.

**Figure 8.** Transverse profiles of maximum velocity in a cross-section of Gage 6 (*Bg* = 0.07 m).

#### *3.3. Effects of Forest on Tsunami Run-up at Actual Scale*

#### 3.3.1. Coastal Topography and Forest Conditions

Coastal forests are widely considered to mitigate tsunami damage on coastal beaches in some countries and regions. In this case study, a numerical simulation was performed to investigate the run-up of tsunami and the effects of forest vegetation on tsunami mitigation. A coastal topography adopted for tsunami simulation is presented in Figure 9 [25,31]. The bed profile of the computed domain consisted of four slopes: S = 1/10, 1/100, 1/50, and 1/500, while the offshore water depth with a horizontal bed was 102 m below the still water level, all parameters for a tsunami wave simulation can be seen in references [25,31].

*P. odoratissimus* is a representative tree of South and Southeast Asia with a complex aerial root structure growing in coastal beach. In addition, *C. equisetifolia* is another representative species that grows densely on sandy beaches. These two species are discussed in evaluating the resistance of vegetation on tsunami propagation. The tree heights *hv* of *P. odoratissimus* and *C. equisetifolia* are considered as 8 m and 11 m, respectively; the reference diameters *bref* of *P. odoratissimus* and *C. equisetifolia* are set as 0.2 m and 0.15 m, respectively; and the density values *N* of *P. odoratissimus* and *C. equisetifolia* are 0.22 trees/m<sup>2</sup> and 0.4 trees/m2, respectively. The depth-averaged equivalent drag coefficient *CD(h)* are variable with vegetation parameters and inundation depth *h*, and they are modified by Tanaka et al. and Thuy et al. [22,26]. The coastal forest starts at a slope of 1/500 (see Figure 9).

In this numerical simulation, the computational domain is discretized into 38,010 triangular meshes, in which the coarsest grid spacing at the sea-side is 20 m, and the finest grid spacing is 2.5 m at the nearshore area. The Manning's roughness coefficient is set as 0.02 for the computational domain. By comparing the maximum water surface elevations in Figure 10, there is good agreement between the simulated and Thuy's results [31].

In non-vegetated areas, as the wave crest approaches a 1/50 slope, its energy compresses and the maximum water surface elevations increase slightly because of shoaling. When the wave crest approaches a slope of 1/500, the maximum water surface elevations decrease due to the accelerated flow velocity on a mild slope. For fully vegetated areas, the maximum water surface elevations increase significantly at the front of the forest due to the reflection of the waves caused by the vegetation. In the forest zone, the maximum water surface elevations decrease drastically. The run-up height above the still water level is 5.10 m with full vegetation, which is lower than the 6.86 m for non-vegetation areas. The results indicate that *P. odoratissimus* effectively protects coastal communities from the destruction by a tsunami wave by reducing wave energy.

**Figure 9.** Schematic topography of a tsunami flume with vegetation effects.

**Figure 10.** Maximum water surface elevations of the first wave without vegetation and with *P. odoratissimus* vegetation [26].

#### 3.3.2. Effects of Forest with a Straight Open Gap on Tsunami Run-up

Figure 11 shows four types of gap arrangements (Cases 1, 2, 3, and 4) in *P. odoratissimus* forest with the simulation conditions of LF = 200 m, BF = 200 m, and BG = 15 m. In this section, tsunami run-up on the sloping beach is considered using a straight gap (Case 1) to investigate the effect on tsunami run-up.

**Figure 11.** Sketch of gap arrangements, (**a**) Case 1, (**b**) Case 2, (**c**) Case 3, (**d**) Case 4.

Figure 12a–c presents the predicted flow structure in Case 1 at different times (600, 660 and 700 s). The flow velocity of the first wave is fast in an open gap zone, while the wave moves slowly in the forest domain due to the resistance from vegetation. The first wave front quickly passes the open gap and spreads out from the exit, and the flow behind the forest gradually becomes uniform, with the exception of the areas close to the gap exit. Figure 13 shows the temporal variation of water surface elevations and flow velocities at four predicted stations (G1, G2, G3, and G4) for Case 1. In this case, the maximum water surface elevations in the presence of vegetation are higher at G1 (the middle of the open gap) and G2 (the middle before the vegetation) than those in the absence of vegetation due the reflection of tsunami waves caused by *P. odoratissimus* forest.

**Figure 12.** Flow structure in Case 1 at different times ((**a**) T = 600 s, (**b**) T = 660 s, (**c**) T = 700 s).

**Figure 13.** *Cont*.

**Figure 13.** Temporal variation of water surface elevations and flow velocities at four predicted stations of Case 1, (**a**) G1, (**b**) G2, (**c**) G3, (**d**) G4.

As the flow passes the forest, the maximum water surface elevations drastically decrease. The maximum water surface elevation values at G3 (the middle behind the forest) and G4 (the middle at the end of the open gap) in the presence of vegetation with an open gap decrease to 80.72% and 82.52% of that for the case of non-vegetation, respectively. The time delays in a tsunami arriving at G3 and G4 are about 45 and 14 s, respectively, when the wave passes through forest with an open gap compared to that of non-vegetated case. The *P. odoratissimus* forests slightly impact the peak velocities at G1 and G2, which reach 4.36 and 4.35 m/s, respectively. Due the resistance effects of *P. odoratissimus*, the flow velocities at G1 and G2 decrease dramatically after reaching their peak. The peak velocity at G3 is decreased by 35.13% and the peak velocity at G4 is increased by 44.04%. The open gap in forest vegetation amplifies negative effects on tsunami propagation even more than in the non-vegetated case.

#### 3.3.3. Effects of Arrangements of Vegetation Patches and Open Gaps

Four types of gap structures are discussed to investigate their potential to increase the energy attenuation of tsunami waves on a vegetated beach. Figure 14 illustrates the structures of a flow field at 660 s for three cases with different gap arrangements, and shows that the flow velocities are larger in open gaps than those in vegetated zones, and the front of the first wave quickly passes the open gap, and the maximum peak flow velocities appear around the gap exit. The furthest distance that tsunami waves reach is variable, but there are differences in flow where structures occur. As presented in Figure 15, the effects of the gap arrangements on the water surface elevations at G3 and G4 are insignificant for four cases, while the different gap arrangements have important effects on flow velocities of tsunami waves.

Compared to Case 1, the peak flow velocities at the G3 station are increased by about 0.15 m/s (4.92%) and 0.21 m/s (6.71%) for Case 2 and Case 4, respectively, and the peak flow velocity is attenuated slightly by about 0.05 m/s (−1.46%) for Case 3. Compared to Case 1, the peak flow velocities at G4 are decreased by about 0.47 m/s (−6.77%), 0.34 m/s (−4.87%) and 0.89 m/s (−12.88%) for Case 2, 3, and 4, respectively. Figure 16 shows the comparisons of the run-up heights of the first tsunami wave for different cases including fully vegetated, non-vegetated, and four types of vegetation arrangements with different gaps. The results indicate that the slight differences of run-up heights occur for different gap arrangements, however, the presence of vegetation significantly attenuates the energy of tsunami waves compared to non-vegetated case. Therefore, it is reasonable to reduce the disadvantages from open gaps by changing gap arrangements.

**Figure 14.** Flow structures at 660 s: (**a**) Case 2, (**b**) Case 3, and (**c**) Case 4.

**Figure 15.** Temporal variation of surface water elevations and current velocities at (**a**) = G3 and (**b**) = G4.

**Figure 16.** The comparison of maximum run-up heights.

3.3.4. Effects of Forest Parameters on Tsunami Run-up

To investigate the effects of coastal forest parameters on tsunami run-up, a sensitivity analysis of different tree configurations, densities, and growth stages was examined in the gap arrangement of Case 4. In this section, Collocation 1 (pure *P. odoratissimus*), Collocation 2 (pure *C. equisetifolia*), Collocation 3 (*C. equisetifolia* as the back vegetation layer of *P. odoratissimus*), and Collocation 4 (*C. equisetifolia* as the front vegetation layer of *P. odoratissimus*) are considered to analyze effects of tree configuration on tsunami propagation. Figure 17 shows the variation of normalized run-up heights (*η*/*η*0) and normalized peak flow velocities (*V*/*V*0) at simulated gauges (G3 and G4) with different tree arrangements, respectively, where *η*<sup>0</sup> and *V*<sup>0</sup> denote the run-up height and the peak flow velocity in the absence of vegetation. The results indicate that *C. equisetifolia* as a new vegetation

layer at the back of existing *P. odoratissimus* (Collocation 3) can minimize the amplification of the peak flow velocity through the open gap, while generating a slightly larger peak flow velocity at the G3 station. The run-up height is higher than that of Collocation 1 (pure *P. odoratissimus*) and smaller than that of Collocation 2. This is because the peak flow velocity is reduced significantly due greater resistance of *P. odoratissimus* when it passes the *P. odoratissimus* vegetation layer, and in back zone of *C. equisetifolia* vegetation with smaller drag force, the velocity increment in open gap is smaller compared to that of other tree configurations.In Collocation 1 (pure *P. odoratissimus*), Figure 18 indicates that normalized run-up height and normalized peak flow velocity at G3 decreased monotonously with increased densities, while the normalized peak flow velocity of G4 first increased (vegetation density ≤ 0.15 tree/m2) and then presented adverse variation. Figure <sup>19</sup> indicates that the increased vegetation height and diameter at different growth stages of *P. odoratissimus* generate positively correlated effects on normalized peak flow velocity at G4, and significant attenuation effects on normalized run-up height and normalized peak flow velocity at G3. Generally, changing forest parameters (e.g. vegetation collocations, densities and growth stages) could be a more cost-effective way to mitigate tsunami damage.

**Figure 17.** Variation of normalized run-up heights and normalized peak flow velocities at G3 and G4 with different tree collocations, (**a**) The normalized run-up heights, (**b**) The normalized peak flow velocities.

**Figure 18.** Variation of normalized run-up heights and normalized peak flow velocities at G3 and G4 with different densities of pure *P. odoratissimus.* (**a**) The normalized run-up heights, (**b**) The normalized peak flow velocities.

**Figure 19.** Variation of normalized run-up heights and normalized peak flow velocities at G3 and G4 with different heights of pure *P. odoratissimus*; (**a**) The normalized run-up heights, (**b**) The normalized peak flow velocities.

#### **4. Conclusions**

Based on the shallow water equations of mass and momentum, an explicit depth-averaged 2D numerical model was developed to simulate wave attenuation and tsunami wave mitigation in a vegetated coastal region. The finite volume method was used to keep conservation of mass in this model. The proposed model was well verified for both breaking and non-breaking solitary wave propagations on a bare sloping beach, and long periodic wave propagation on a partially-vegetated sloping beach, as shown by comparison of experimental data. Overall, the present numerical model has been shown to accurately account for wave propagation and transformation in vegetated and non-vegetated sloping beaches.

In the case of full *P. odoratissimus* forest on an actual-scale beach, the water surface elevations at the front of the forest belt increased due to the effect of reflected waves caused by vegetation. Here, *P. odoratissimus* forest played a positive role in reducing maximum water surface elevations through the vegetated zone, and as such the maximum wave run-up height of the tsunami was effectively reduced. When a forest has an open gap, the negative effects can be reduced by choosing a beneficial arrangement of vegetation and open gaps. Sensitivity analysis was also carried out by changing tree collocation, tree densities, and growth stage, with results showing that they had important impacts in reducing the peak flow velocity of the tsunami wave at the open gap exit.

Overall, these study findings provide cost-effective natural strategies to improve the effectiveness of *P. odoratissimus* forests with an open gap by changing vegetation arrangements and vegetation parameters (collocations, tree densities, and growth stage), which are important in the design and establishment of vegetated bioshields against tsunami hazards.

**Author Contributions:** H.Z. contributed to the establishment of model, significantly to analysis and manuscript preparation. M.Z. contributed to the conception of the study and the final editing. T.X. and J.T. contributed to the collection of experimental data.

**Funding:** This work was supported by the National Nature Science Foundation of China (51779039, 51879028), the Wetland Degradation and Ecological Restoration Program of Panjin Pink Beach (PHL-XZ-2017013-002), the Fund of Liaoning Marine Fishery Department (201725).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Water* Editorial Office E-mail: water@mdpi.com www.mdpi.com/journal/water

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18