**A Review on Hydrodynamics of Free Surface Flows in Emergent Vegetated Channels**

**Soumen Maji 1, Prashanth Reddy Hanmaiahgari 2,\*, Ram Balachandar 3, Jaan H. Pu 4,\*, Ana M. Ricardo <sup>5</sup> and Rui M.L. Ferreira <sup>6</sup>**


Received: 27 February 2020; Accepted: 16 April 2020; Published: 24 April 2020

**Abstract:** This review paper addresses the structure of the mean flow and key turbulence quantities in free-surface flows with emergent vegetation. Emergent vegetation in open channel flow affects turbulence, flow patterns, flow resistance, sediment transport, and morphological changes. The last 15 years have witnessed significant advances in field, laboratory, and numerical investigations of turbulent flows within reaches of different types of emergent vegetation, such as rigid stems, flexible stems, with foliage or without foliage, and combinations of these. The influence of stem diameter, volume fraction, frontal area of stems, staggered and non-staggered arrangements of stems, and arrangement of stems in patches on mean flow and turbulence has been quantified in different research contexts using different instrumentation and numerical strategies. In this paper, a summary of key findings on emergent vegetation flows is offered, with particular emphasis on: (1) vertical structure of flow field, (2) velocity distribution, 2nd order moments, and distribution of turbulent kinetic energy (TKE) in horizontal plane, (3) horizontal structures which includes wake and shear flows and, (4) drag effect of emergent vegetation on the flow. It can be concluded that the drag coefficient of an emergent vegetation patch is proportional to the solid volume fraction and average drag of an individual vegetation stem is a linear function of the stem Reynolds number. The distribution of TKE in a horizontal plane demonstrates that the production of TKE is mostly associated with vortex shedding from individual stems. Production and dissipation of TKE are not in equilibrium, resulting in strong fluxes of TKE directed outward the near wake of each stem. In addition to Kelvin–Helmholtz and von Kármán vortices, the ejections and sweeps have profound influence on sediment dynamics in the emergent vegetated flows.

**Keywords:** turbulence; emergent vegetation; flexible vegetation; rigid vegetation; coherent structures; shear layer

#### **1. Introduction**

Vegetation is ubiquitous in rivers, estuaries, lake shores and some coastal areas [1]. It can be of different types including submerged, floating and emergent, whose incidence in a channel is schematized in Figure 1. Plants having varied stiffness, density, flexibility and height influence the water flow in different ways [2,3]. Its spatial distribution is seldom uniform [4]. Vegetation is frequently arranged in patches [5,6]. The presence of vegetation in water bodies has major impacts on flow hydrodynamics, including mean flow and variables describing turbulence. Certain important features like flow rate, changes in the bed and sediment carrying capacity of the stream are influenced by riparian plantation and its interactions with hydrodynamics.

**Figure 1.** Different types of vegetation in open channel flow.

Parameters such as increased length, density and height of vegetation result in increased roughness, increased water levels and decreased velocities. In some cases vegetation offers protection against bank erosion, tsunamis and high waves. Experimental results show that the presence of vegetation leads to the generation of turbulence and additional drag forces. It is well known that the vertical distribution of flow velocity depends on the density of vegetation. The increase in stem density decreases the flow velocity through the interior of vegetation and increases the velocity in the regions without vegetation. Lichtenstein [7] classified aquatic plants into four types: algae, floating-leaved plants, submerged plants and emerged plants as shown in Figure 2a.

(a)

**Figure 2.** *Cont.*

**Figure 2.** (**a**) Classification of aquatic plants; (**b**) emergent vegetation types based on plant stiffness.

Vegetation may be broadly divided into rigid and flexible, as shown in Figure 2b. The effect of flexible vegetation is quite different from that of rigid vegetation. The rigid vegetation causes the flow to produce separation, wakes and eddies which in turn dissipate the flow energy [8]. Figure 2 shows definition of emergent, rigid and flexible vegetations, in which, *HV* = height of rigid vegetation, *hc* = height of bent vegetation and *H* = water depth.

Given the complexity of emergent vegetated flows with ecological, geomorphologic and hydrodynamic effects, a great deal of research work including sophisticated field, laboratory and numerical studies have been undertaken to understand it. Significant advances have been achieved in the field of emergent vegetation hydrodynamics with respect to horizontal and vertical structures, [3,4,8–13], drag and frictional characteristics, [9,10,14–19]. Furthermore, with increasing computer power and improvements in numerical modelling, impressive advances were made in the study of the origin and development of coherent structures in the emergent vegetated flows [11,12]. These numerous investigations require a consolidated review to help researchers working on the hydrodynamics of emergent vegetated flows for deeper understanding. Therefore, the aim and scope of the paper is to make a synthesis of the state of the art relevant research works on the hydrodynamics of emergent vegetated flows, with emphasis on mean flow, turbulence and drag. In line with these objectives, present review highlights important research contributions of emergent vegetation hydrodynamics including complex interactions between turbulence and vortex structures.

The review of research findings on open channel flows with emergent vegetation is further divided into: (i) vertical structure of flow field (ii) velocity distribution, 2nd order moments, and distribution of turbulent kinetic energy (TKE) in the horizontal plane, (iii) horizontal features which include wake and shear flows and (iv) drag effect of vegetation on the flow. In addition, future directions are deliberated for researchers working on the emergent vegetated flows, and finally conclusions of the paper are presented at the end.

#### **2. Hydrodynamics of Emergent Vegetated Flows**

Vegetation extending from the bed all the way beyond the free surface is called emergent vegetation. Good examples of emergent vegetation include jute plantations, mangroves, and often partially submerged trunks of palm trees in flood water. Studying emergent vegetation is sought after because it modifies velocity vertical profiles and turbulence intensities. Commonly, emergent vegetation consists of relatively stiff stems, which are usually circular in cross section and without leaves below the water surface [3]. Therefore, emergent rigid cylindrical stems have been extensively employed in research works as a proxy for natural stems.

In unbounded open channels, the mean flow through groups of emergent vegetation is decreased. The emergent canopy has a significant role on coherent structures and accompanying mass and momentum fluxes at the boundary between emergent vegetation and open channel. In particular, emergent canopy patch density in the open channel flow affects mass and momentum transfer, roughness, sedimentation, velocity of flow, bed shear, turbulence quantities, biological processes and aquatic life within the canopy. Additionally, characteristics of flow through emergent vegetation consist of significant velocity gradients (in transverse, stream wise and vertical directions) and drag discontinuity at the interface resulting in the formation of a shear layer between the vegetation stems and the flow exterior of the vegetation. The shear layer in turn generates large coherent vertical structures due to Kelvin–Helmholtz instability. These coherent structures grow downstream due to creation of vortex-street, vortex pairing and eventually dissipate into the flow.

#### *2.1. Vertical Structure*

The vertical structure of flows within the stems of emergent vegetation is divided in three layers: (i) the reach close to the bottom, where the flow is highly 3D due to interaction with the bed; (ii) the layer close to the free surface, which is affected by the oscillations of the free surface; and (iii) the intermediate layer, sufficiently away from the bed and from the free surface, where the flow is controlled by the vertical stems and the flow properties are approximately constant in the vertical direction [9–11,20] as shown in Figure 3.

**Figure 3.** Vertical structure of flow in emergent vegetation.

Stoesser et al. [11] numerically studied three types of spacing between cylinders on the flow field. Very small vertical velocities are observed in the large spacing cases except just behind the cylinder where relatively large values of vertical velocity indicates considerable upward movement of fluid observed just behind the cylinder. However, upward and downward movement of fluid is observed in the case of low spacing between the cylinders. Contours of the time-averaged streamwise velocity and streamlines at about half depth for the three different vegetation densities are investigated. Flow separation occurs at approximately 95◦ and a relatively large recirculation region comprising two counteracting vortices that are about the size of the cylinder diameter are observed in low and medium vegetation densities. In high vegetation density, flow separates considerably later and the recirculation region behind the cylinder is much smaller.

Ricardo et al. [3] addressed the influence of the longitudinal variation of the stem areal number density on the spatial distribution of the turbulent flow variables at the inter-stem scale. They concluded that Reynolds stresses are not sensitive to local spatial gradients of the stem distribution, being determined by the local number of stems per unit area. However, form-induced stresses depend on the local number of stems and its spatial distribution. Ricardo et al. [4] identified a spatial memory of the flow manifested by the upstream generated complexity of the flow that subsists further downstream and is combined with the complexity locally generated, resulting in the increased spatial variability relative to that of a flow with a uniform stem distribution.

Chang and Constantinescu [12] numerically investigated the flow and turbulence structure interior and exterior of a circular array of emergent rigid vegetation patch by varying the solid volume fractions and relative cylinder diameter arrays. These studies were performed without considering a different bedforms effect, which could be a crucial factor in flows through natural bedforms as studied by Pu et al. [21]. The flow was found to be highly three dimensional in the array when solid volume fraction (SVF) was more than 0.1 and this effect was relatively small when the magnitude of SVF was less than 0.025. Von Kármán vortex streets were developed for cylinder Reynolds number greater than 10,000. Furthermore, Chang and Constantinescu [12] found that at sufficiently high SVFs such as 0.2, von Kármán wake cylinder interactions become important resulting in the formation of von Kármán vortex street. In low SVFs such as 0.023, wake and cylinder interactions were weak and flow advects through the cylinders without a vortex street. However, shear layers are observed irrespective of formation of von Kármán vortex streets. For larger values of SVF, cylinders in the downstream imped the formation of strong vortices.

Liu and Shan [22] used an analytical model to predict longitudinal profile of depth averaged streamwise velocity for simulated emergent cylindrical vegetation and used an experimental setup to validate the proposed model for short and long arrays of vegetation. Tong et al. [23] investigated the effect of emergent rigid vegetation in a *y* shaped confluence channel and concluded that non submerged vegetation is responsible for changing the flow structure. Due to the vegetation, middle of the non-vegetated area in the vicinity of tributaries showed existence of high velocity region, and the disappearance of the secondary flow. The flow in non-vegetated areas was more intense and turbulent kinetic energy of the non-vegetated area was lesser than the vegetated area.

#### *2.2. Mean Flow and Distribution of Turbulent Kinetic Energy (TKE) in the Horizontal Plan*

Ricardo et al. [8] carried out a laboratory work simulating rigid and emergent cylinders, with longitudinally varying density where 2D instantaneous velocity maps were measured with a Particle Image Velocimetry (PIV) aiming at characterizing the mean flow field at the inter-cylinder and patch scales. Figure 4a,b exemplifies a horizontal map of longitudinal velocity and out-of-plane vorticity map for a patch of 1200 cylinders/m<sup>2</sup> in a 3.5 m-long reach populated with rigid cylinders and with longitudinally varying density of cylinders.

The mean flow field within vegetation covered areas is heterogeneous at large scales, as reach or patch scales, with zones of low longitudinal velocity at the vegetation elements wake alternated by high longitudinal velocity zones between those elements [8,24]. These low/high velocity patterns are observed independently of the density and distribution of the vegetation elements.

**Figure 4.** (**a**) Time-averaged longitudinal velocity map (m/s); (**b**) out-of-plane vorticity map (m<sup>−</sup>1); (**c**) rate of production of TKE (m<sup>2</sup> s<sup>−</sup>3), at longitudinal position P8 and at 3.1 cm above the bed corresponding to 60% of the flow depth. (Source: Ricardo et al., 2014).

Relative to the vorticity distribution, a repeating pattern of paired vortexes caused by the unsteady separation of the flow on the cylinders is observed. These quasi-symmetric high vorticity patterns behind the cylinders identify von Kármán vortex streets. Comparing vorticity maps for the different density of cylinders, the conclusion is that the cylinders induce a regular structure of vortex patterns independently space between cylinders. The main difference is that the space necessary to fully develop the vortex pattern is strongly reduced at smaller distances between cylinders. At higher cylinder densities, the vortex pattern is forced to compress due to the proximity of the next cylinder, while at positions with lower cylinder density cases the vortex pattern has space to develop a vortex street.

Ricardo et al. [8] performed laboratory experiments on rigid emergent vegetation aimed at the characterization of the key terms of the TKE equation. Three dimensional laser doppler anemometry (LDA) and two dimensional particle image velocimetry (PIV) measurements have been employed. They characterized the spatial distribution of turbulent production, convective rate of change of TKE and turbulent diffusion and the dissipation rate of TKE was also computed. They observed that the production of TKE is mostly associated with vortex shedding from individual stems (Figure 4c). The magnitude of the rate of production is higher in the wake region and negative production of TKE was identified between close cylinders, associated with strong local accelerations in the flow field. Turbulent transport is particularly important along the von Kármán vortex street and convective rate of change of TKE and pressure diffusion are most relevant in the vicinity of the cylinders. The rate

of dissipation was found to increase with the stem areal number density. Ricardo et al. [8] observed that the rates of production and dissipation are not in equilibrium in the inter-stem space, revealing important interactions of turbulence with mean flow and pressure field and turbulent transport of TKE. The cumulative effect of convection and turbulent transport of TKE is the generation of background turbulence, i.e., non-locally generated turbulence.

Ricardo et al. [4] investigated the impact of the background turbulence caused by randomly placed cylinders on the vortex shedding regime. The strong background turbulence, generated by the randomly placed cylinders in the array, causes a premature vortex detachment increasing the shedding frequency. It was concluded that the background turbulence acts on the faster loss of vortex coherence as they travel downstream. Furthermore, Ricardo et al. [4] observed that the impact of the background turbulence on decay of the longitudinal vorticity flux is small, the latter is mainly caused by the vorticity cancellation imposed by the local distribution of neighbouring cylinders.

#### *2.3. Horizontal Structure*

Density, shape and arrangement of the stems in the emergent vegetation patch have a dominant role on the mean flow and turbulence structures in the interior and exterior of the patch in horizontal direction.

The mean flow is largely dominated by the wakes which tend to reduce the flow between the bodies and this effect is larger than the inviscid kinematic effect of the bodies blocking the flow. The difference between Lagrangian and Eularian velocities with different void ratios is one of the aspects of horizontal structure. At the upstream of the array of cylinders, the flow is irrotational and slows down due to the blocking effect of the entire patch and the drag force from the cylinders. The flow bleeds through the vegetation with low velocity and high velocity flow takes through the region without vegetation. Shear layers form on either sides of the vegetation patch as shown in Figure 5. At what distances upstream and downstream of the flow is influenced by the array of cylinders is a function of vegetation density, width of the patch, Reynolds and Froude numbers.

Stoesser et al. [11] carried out large eddy simulations (LES) of flow through the canopies with solid volume fractions 0.016, 0.063 and 0.251. In their analysis, the streamwise turbulence intensity profiles were significantly different from unobstructed channel flow not only in shape but also in magnitude. Higher streamwise turbulence intensities were observed in high vegetation density cases. Similarly, normalized vertical turbulence intensities have higher magnitudes for higher density vegetation. The turbulence in the high vegetation density case is rather generated by high streamwise velocity gradients, while von Kármán vortex shedding occurs in low density cases [11].

Nicolle and Eames [25] studied the effects of the patch of emergent and rigid cylinders on the flow field in a uniform air flow using direct numerical simulations (DNS). Nicolle and Eames [25] found that the von Kármán vortex-street behind the circular cylinder array is similar to that of flow past a single cylinder. The most obvious limitation of their study is the implementation of 2D model to study the 3D flow field where additional drag is caused by vortex stretching.

Nepf [26] concluded that emergent vegetation develops two different types of turbulent scales i.e., canopy scale turbulence and stem scale turbulence in which canopy scale vortices control the high level of turbulent diffusion in the exchange zone. In general, the turbulence in the interior of emergent vegetation is less than that of the outside the vegetation in free surface flows. However, the decrease in turbulence due to decreased mean velocity near the bed is too small as compared to the additional turbulence generated by stem vortices. White and Nepf [13] and White and Nepf [27] considered horizontal flow structure and studied the momentum transfer at the interface between the flow inside the emergent vegetation and the region outside the vegetation. Shear layer is created due to the drag discontinuity at the interface and coherent vortices are generated in the shear layer due to Kelvin–Helmholtz instability. The exchange of fluxes between the open channel flow through the vegetation and the region outside the vegetation is dominated by the energetic shear layer vortices. Vortices in the shear layer are found to be continuously growing in the downstream direction

predominantly through vortex pairing. In addition a quasi-coherent vortical structure with two times the spacing of the shear layer vortices is generated by the interaction of paired shear layer vortices.

)ORZGLUHFWLRQ

**Figure 5.** Horizontal structure of the flow through a vegetation patch.

Meftah and Mossa [28] studied channel partially obstructed by emergent vegetation with an emphasis to understand the effect of contraction ratio. They observed the formation of a shear layer immediately next to the interface, followed by an adjacent free stream region of full velocity flow. In addition, the transversal flow velocity profile at the obstructed and unobstructed interface was analyzed. The experimental results demonstrated that the contraction ratio significantly affects the flow hydrodynamic structure. Finally, a general modified log-law was proposed to describe the typical transversal profile of the mean velocity and its application was in good agreement with analytical solutions.

Naot et al. [29] used computational fluid dynamics (CFD) with the *k* − model to simulate hydrodynamic behaviour of compound channel with rigid and emergent vegetation on the flood plain, they found that flow behaviour depends on shading factor and reference length of the vegetation. A two dimensional LES model was used by Ikeda et al. [30] to get a relation between the square of flow velocity and the waving motion of a flexible emergent (where plant height is equal to water depth) plant. The interesting observation in this study is the motion of the plant increases the frequency of vortices whereas increased strength of vortices increases the motion of plants.

Stoesser et al. [11], and Kim and Stoesser [31] concluded that von Kármán-type vortices, Kelvin–Helmholtz instability and vortex shedding are clearly visible in sparse vegetation. As the vegetation density increases, the influence of vortices that are shed from upstream cylinders affects more irregular shedding behavior of downstream cylinders. Increase in vegetation patch density increases the vorticity magnitude which is accompanied by the linear increase of vertical vorticity. Furthermore, the frequency of vortex shedding increases as the patch density increase.

Chang and Constantinescu [12] numerically found that coherent structures formed in the separated shear layer (SSL) are larger with SVF = 0.2 as compared to SVF = 1.0. For low SVF values, the separated shear layer (SSL) was found to be longer and eddies in the SSL have more energy as compared to that of the SSL in high SVF patch. Furthermore, the behaviour of eddies in SSL of low SVF patch resembles that of a mixing layer. Von Kármán Street was formed only in flows with SVF values higher than 0.2 and in which vortical structures were impeded and, therefore, not producing the regular wake vortices. For SVF values less than 0.05, shedding of wake vortices behind a cylinder was strongly disturbed by the vortices shedding from other cylinders in the vicinity. Chang and Constantinescu [12] concluded that for SVF = 0.023, shedding behind an individual cylinder in the patch is similar to that of an isolated cylinder in the flow.

Anjum and Tanaka [32] used discontinuous and vertically doubled layered patches and simulated it numerically for varying vegetable density and patches. Large-scale turbulence followed by saw-tooth distribution within the patches and low turbulence in the non-vegetation regions was observed. Maximum of 13% turbulence intensity for dense vegetation arrangements and maximum of 9% turbulence intensity for sparse vegetation was observed. They concluded the presence of non-uniform flow through discontinuous and double layered vegetation patches. Yamasaki et al. [33] studied the hydrodynamics of flow-vegetation interactions using a computational fluid dynamics model.

#### *2.4. Drag and Frictional Characteristics*

In open channel flows without vegetation, streamwise velocity distribution along the depth is a function of Reynolds shear stress, whereas velocity vertical profile in the vegetated channel is a function of frictional drag imposed by the vegetation since the vegetation drag is far greater than the frictional drag at the bed in high-density vegetated flows [10,34]. The drag force exerted by the flow on emergent vegetation can be estimated based on the available knowledge of flow around a circular cylinder.

Many studies are available in the literature to predict or measure frictional resistance in emergent vegetation conditions [9–11,14–16,35,36]. To compute the flow resistance in a vegetated channel, the drag caused by a single stem and the addition of individual drags of all the stems in that canopy, the blockage effect of the canopy in addition to the already existing bed roughness must be considered. Fathi-Maghadam and Kouwen [16] were among the first researchers who studied emergent flexible vegetation effect on the drag coefficient and in-turn onManning's roughness coefficient. Fathi-Maghadam and Kouwen [16] experimentally analysed the effect of proposed non-dimensional variables on *Cd* and concluded that (i) Manning's roughness proportionally increases with an increase in the square root of flow depth and increase in density of vegetation, and (ii) Manning's roughness inversely proportional to the average flow velocity and vegetation flexibility.

In the riverine environment, however, canopies occur in finite patches and make a random angle with the direction of flow; for example, patches of mangroves in Purna River, India and Daintree River Queensland, Australia as shown in Figure 6. In these scenarios, flow is divided into flow passing through interior and exterior of the vegetation. According to Mitul et al. [18], mangroves are the densest canopies in coastal water with stem diameters between 4 cm to 9 cm and their solid volume fraction is around 0.55.

**Figure 6.** Patches of vegetation (highlighted in circles) present in the river systems (Background image is sourced from Google Earth).

Maji et.al. [37] investigated the turbulence characteristics in open channel flow with an emergent, sparse and rigid vegetation patch. They observed that the nature of time-averaged velocities along streamwise, lateral and vertical directions are not consistent upto half the length of the patch from its leading edge; however, velocity profile becomes uniform after that length and their results are similar to various un-obstructed uniform and non-uniform flows studied by Pu et al. [38]. Inside the patch, the magnitude of vertical normal stress is small as compared to the streamwise and lateral normal stresses and Reynolds shear stress decreases along the centreline of the patch in the downstream. Furthermore, flow velocity decreases within the vegetation, which can lead to deposition of sediment which accelerates with the time and acts like a catalyst to the larger morphological changes in the river [37,39].

Heidari et al. [40] studied the drag created by a slender emergent cylinder in open channel flow and proposed that wake stability parameter significantly affects the drag on the cylinder away from the bed. Kothyari et al. [41] developed a drag model applicable to emergent plants in which drag coefficient decreases as the stem Reynolds number increases. Luna et al. [42] defined a term called 'global flow resistance' similar to Chezy's coefficient, however the global flow resistance is a combination of bed shear stress and vegetation drag on the flow. Luna et al. [42] assessed various global flow resistance models for emergent vegetation and found that Baptist [43] model accomplishes the bestfit. Takemura and Tanaka [44] investigated the fluid structures in a colony type of emergent vegetation and coefficient of drag was experimentally measured. Their results restated that vortex patterns, strouhal number and drag coefficient of colony are function of *L*/*D* where *L* is cross stream spacing between neighboring cylinders and *D* is cylinder diameter. Takemura and Tanaka [44] termed a von Kármán vortex-street behind an individual cylinder as primitive von Kármán vortex street (PKV) and its instability governs large scale vortical structures.

Tanino and Nepf [9] examined the canopy drag as a function of*Rep* = (*UD*/ν, where *U* = double averaged pore velocity, υ is kinematic viscosity and *D* is stem diameter) which is in the range of 25–685 and solid volume fraction ranging from 0.091 to 0.35. They found that the *Cd* value of a canopy patch is inversely proportional to *Rep* and proportional to solid volume fraction. Average drag per unit cylinder length non-dimensionalized by the product of *U* and dynamic viscosity (μ) is found to be a linear function of *Rep*. This linear relation was further investigated by Ferreira et al. [10], who proposed that the increase of the solid fraction generates a momentum sink observed as an increase in the gradients of dispersive stresses. This may justify the increase of the drag coefficient with the increase of solid fraction at intermediate values of the latter. Liu and Zeng [45], however, suggest that, in general, the drag coefficient decreases with the solid fraction. Another important observation is that the drag of an individual cylinder is not a function of solid volume fraction. Chen et al. [17] studied the wake region behind a circular vegetation patch and found that for low flow blockages sediment deposition occurs, whereas for high flow blockages erosion takes place in the wake region.

Musleh and Cruise [36] combined the effects of flow depth, velocity, stem diameter, lateral spacing and longitudinal spacing into a non-dimensional variable called relative density ratio (ratio of total surface area of the vegetation to bed cross-sectional area) that is linearly related to Darcy's friction factor and mannings coefficient. Musleh and Cruise [36] explained that for rigid vegetation, flow resistance of vegetation non-linearly increases with increasing flow depth especially when vegetation density is high, whereas flow resistance found to be nonlinearly decreasing with increasing mean flow velocity for flexible emergent vegetation flows. In addition, they found that the arrangement of rigid emergent cylinders plays an important role on the flow resistance and friction factor linearly varies with the relative density ratio. Musleh and Cruise [36] found that lateral spacing and stem diameter of cylinders have a more significant effect on the flow resistance as compared to the streamwise spacing of cylinders for a given stem density. The important observation in this experimental study is the drag effect of stem diameter is more than the transverse and streamwise spacings of cylinders. Furthermore, Musleh and Cruise [36] found that lateral spacing between the rods significantly affects the drag on rods. Wu et al. [19] experimentally found that the drag coefficient of emergent rigid vegetation increases with the increasing flow Reynolds number.

Chen et al. [14] used a tapered emergent flexible vegetation plant model to derive an analytical solution for flow resistance, bending behaviour and potential rupture locations. In this method, power law was used to determine the velocity distribution which is responsible for drag force on the rigid as well as flexible stems. Analytically computed deflections were validated with experimentally measured deflections. Finally, Chen et al. [14] concluded that the location of rupture depends on the tapering angle and the base diameter of the stem. Their results were consistent with earlier research that drag is inversely proportional to the average flow velocity.

Stoesser et al. [11] found that value of *Cd* is linearly proportional to canopy density. In highly dense canopies, the downstream cylinders experience a drag force which is slightly smaller than for an isolated cylinder, but the lift force is much larger than the isolated cylinder because downstream cylinders are sitting in the wakes of two upstream cylinders [25]. In sparse canopies (SVF 1.0) upstream cylinders experience a drag force similar to that of an isolated cylinder but the lift force, while still relatively small is determined by the position of the cylinders in downstream. Cylinders located in the wake of other cylinders experience larger drag and lift forces as a result of the accelerated flow induced between the upstream cylinders. For intermediate void fractions, the force on the individual cylinders is steady and the lift force acting on the whole patch is negligible. The aggregate drag force on the canopy is dominated by the contributions from the upstream rows of cylinders, which experience a drag force comparable to that of a single cylinder. For high void fractions, the lift forces are significant and are mainly induced by the wake of the whole array. The effect of increasing void fraction leads to a slight increase in the scatter of the forces on the upstream cylinders and significantly increases the magnitude of the wake forces on the downstream cylinders.

Stoesser et al. [11] also studied the forces on cylindrical stems in a canopy and derived the following conclusions. Pressure drag, friction drag components are high and bed shear is low on the single cylinder in high patch density as compared to low patch density. In general the drag forces, pressure drag constitutes about 90% and remaining forces are only 10% (frictional drag is around 7% and bed shear is around 3%). Bed shear stress decreases with increasing patch density. The friction drag is almost constant regardless of vegetation density. Chang and Constantinescu [12] found that for numerical simulations, drag coefficient of a cylinder in the array is a very important parameter. They have also added that combined drag parameter for porous cylinders increases as frontal area per unit volume of patch increases.

Rooijen et al. [46] used an isolated emergent cylinder to assess canopy drag and obtained good agreement in drag computation by their proposed analytical model. In the case of a submerged cylinder, the authors suggested that their results are consistent with isolated cylinder theory. Wang et al. [47] used an experimental model to investigate the drag coefficient due to rainfall over vegetation. Experiments were conducted with varying stem density and applied water to obtain a link between friction slope and local kinetic energy head induced by steady non uniform flow. Spatial variation of drag coefficient exhibited a monotonic reduction during rain or non-monotonic hump formation when rain was absent.

Shan et al. [48] examined the effect on the change of drag force and velocity inside mangrove forest model for different vegetation stem arrangements. They found that the force at root is most dominant, randomly distributed branches inside a tree may experience more or less same force, and however, drag force increases with tree density. Razmi et al. [49] used a field scale large eddy simulation model by using distributed-drag approach and velocity-dependent models to study the effect of canopy reconfiguration on canopy posture and canopy drag. This model removed the flaws in the previous research of not capturing vertical spread in peak Reynold stress that is associated with the movement of the canopy interface.

#### *2.5. Vegetation in Bank Protection and Sediment Transport*

The effect of riparian vegetation on river morphology is well documented [50,51]. The role of vegetation density on improved riverbank stability is an emerging topic and important for further research [52]. River bankside vegetation can effectively improve bank strength and stability by mechanical, hydrologic, and hydraulic effects [53–55] and hence it may be considered essential to plant vegetation at erosion-prone locations on the riverbank [56]. Mechanical method comprised of root reinforcement along with flow resistance increment [53,57–59]. Hydrologic method protects the river bank from current induced erosion by reduction of soil moisture and pore water pressure, and increasing the root-induced compaction [60]. By contrast, the hydraulic method causes increase in hydraulic roughness and toe erosion reduction [61–63]. Interconnection between soil and root by apparent cohesion can increase bank stability [64]. But plant root can also negatively affect stability if extra weight of the vegetation increases downward force causing bank failure [58,65]. In addition, vegetal cover can reduce the danger of river course change and braid formation, which helps in maintaining river morphology [59,66–70]. Yu et al. [71] used the bank stability and toe erosion model (BSTEM) to find the effects of different root conditions on channel bank strength and verified that roots are good at reinforcing the unconsolidated banks and erosion control, hence proving the effectiveness of riverbank vegetation.

Yang and Nepf [72] proposed that near bed turbulent kinetic energy is better than bed shear stress for indicating bed load transport in vegetated flows. They proposed a turbulent kinetic energy based Einstein–Brown equation and satisfactorily compared the bed load transport with experimental results. Armanini and Cavedon [73] used a probabilistic/deterministic model to estimate the flow field changes induced by vegetation and associated bed load transport by refining their previous bed load model [74] and recalibrating [75] parameters. Their model shows good resemblance with experimental observations.

#### **3. Further Research Prospects**

Most of the research studies available are based on using rigid cylinders as vegetation stems, but real vegetation has a different shape, roughness and bending stiffness, and therefore laboratory experiments may not substitute the field investigations and this could be the focus in the future. In fact, a few researchers are working towards this achievement. Plant morphological studies were conducted with real vegetation for uniformly vegetated flow [2,76–79]. But there is a lack of dedicated study for partly vegetated flow with real vegetation. Caroppi et al. [80] are carrying out pioneering research

by using naturally grown plants in an open channel for the study of reconfiguration and mixing effect of flow. In the real channels, aquatic seeds grow into aquatic plants and further experimental studies are required to study the transient flow field during aquatic plants growing stage. Similarly, capabilities of numerical simulations must be extended to field plants and scenarios. The feedback between flow field and sediment dynamics in the interior and exterior of the vegetation is beginning to be recognized as vital. Such research now provides an outstanding background for beginning to address areas of greater complexity that will enable a better understanding of these important natural features. In this context, studies conducted by [81–83] on threshold magnitude of velocity fluctuations and duration of velocity fluctuations for particle entrainment from the bed are useful for developing new sediment transport theories based on the energy approach in emergent vegetated flows. It is recognized that sweep events play a dominant role in dislodging bed particles. Flow conditions and the solid volume fraction of emergent vegetation affect the individual contributions of sweep and ejection coherent structures in the interior and wake regions. Therefore, the direct relation between the sediment transport and instantaneous coherent flow structures in vegetated flows are to be probed further. The effect of patch size of aquatic plants has not yet been fully investigated and elaborated clearly. As per the authors' knowledge, the effect of changing the size of the patch in hydrodynamics and sediment transport through vegetation was studied only by [84] and gave an idea of the minimum size of a patch for triggering the flow interaction but it required the study of detailed turbulence characteristics for further investigation.

#### **4. Conclusions**

The additional drag created by vegetation resulting in a decrease in discharge in irrigation channels across the world has prompted researchers to characterize and quantify the effect of emergent vegetation on the flow field. Furthermore, the advent of accurate invasive and non-invasive velocity-measuring devices has accelerated the research in this area since early 2000. It was well accepted that the presence of emergent vegetation in the rivers and open channel flows has a significant role on mean velocity, turbulence levels, coherent structures and the accompanying mass and momentum fluxes at the boundary between emergent vegetation and an open channel. Decreased flow velocity within the vegetation leads to the deposition of sediment which accelerates with the time and acts like a catalyst to the larger morphological changes in the river. Another important and consistent research finding is that universal logarithmic law is not applicable for spatially and time averaged velocity vertical profile within the flow region of emergent vegetation.

Experimental evidence shows that the magnitude of vertical normal stress is small as compared to the streamwise and lateral normal stresses within the vegetation patch. The streamwise and vertical turbulence intensities, vorticity and vortex shedding have higher magnitudes for higher solid volume fraction of the vegetation. Investigation of spatial distribution of the turbulent flow variables at the inter-stem scale demonstrated that Reynolds stresses are not sensitive to local spatial gradients of the stem distribution, whereas form-induced stresses depend on the local number of stems and its spatial distribution. The production of TKE is mostly associated with vortex shedding from individual stems and therefore the magnitude of the rate of TKE production is higher in the wake region. Another important and consistent research finding is that the rate of TKE production and dissipation is found to be increasing with the increasing stem areal density. Few researchers studied the momentum transfer at the interface between the flow region inside and outside the emergent vegetation and observed the formation of coherent vortices in the shear layer due to Kelvin–Helmholtz instability. Researchers found the motion of the flexible emergent plant increases the frequency of vortices whereas the increased strength of vortices increases the motion of plants.

The coefficient of drag of the emergent vegetation patch was found to be inversely proportional to the cylinder Reynolds number, flow velocity and vegetation flexibility. However, the coefficient of drag was found to be proportional to the solid volume fraction, flow depth, and bed slope. It was also established that for low flow blockages, Kelvin–Helmholtz vortices are dominant, whereas for high flow blockages in which the vegetation patch acts like a bluff body von Kármán vortices are dominant. It was established that stem diameter of cylinders has a more significant effect on the flow resistance as compared to the streamwise and lateral spacing of cylinders for a given stem density.

**Author Contributions:** S.M.: writing—original draft preparation, funding acquisition; P.R.H.: writing—review and editing, project administration; R.B.: writing—review and editing, data curation; J.H.P.: writing—review and editing, supervision; A.M.R.: writing—review and editing; R.M.L.F.: writing—review and editing. All authors contributed to the work. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

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

#### **References**


© 2020 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* **Evaluation of Ecosystem Services in the Dongting Lake Wetland**

#### **Li Ma 1, Ruoxiu Sun 1, Ehsan Kazemi 2, Danbo Pang 3, Yi Zhang 4, Qixiang Sun 5, Jinxing Zhou 1,\* and Kebin Zhang 1,\***


Received: 8 October 2019; Accepted: 29 November 2019; Published: 5 December 2019

**Abstract:** The Aeronautical Reconnaissance Coverage Geographic Information System (ArcGIS) 10.2 and Integrated Valuation of Ecosystem Services and Trade-offs (InVEST) model are used to comprehensively evaluate ecosystem services in the Dongting Lake Wetland, focusing on water yield, soil conservation, carbon storage, and snail control and *schistosomiasis* prevention. The spatial and temporal variations of these services, as well as their variations between different land use types in a period of 10 years from 2005 to 2015, are investigated, and the value of such services is then estimated and analyzed. The results of this study show various temporal and spatial trends in the ecosystem services, such as (1) the overall increase of all these services during the study period (although significant in some services, such as *schistosomiasis* patient reduction, by 86.8%; and, very slight in some others such as soil conservation, only by 0.02%); (2) different orders of the services values that are based on different land use types; and, (3) the temporal changes in the proportion of the values of different ecosystem services with respect to the total services value. Besides, it is concluded that the evaluation of ecosystem services of a certain wetland is heavily dependent on the characteristics of the area where the wetland is located, and the assessment indicators and methods should be selected based on such characteristics through the analysis of the results and a comparison with the findings of literature.

**Keywords:** InVEST model; wetland; ecosystem service assessment; value analysis; *schistosomiasis* prevention

#### **1. Introduction**

Ecosystem service is the natural environmental condition and utility that the ecosystem and ecological process form and maintain [1,2]. The millennium ecosystem assessment (MEA) has found that 60% of global ecosystem services have deteriorated, with human activity being one of the main causes of this adverse change [3,4]. Wetlands, as important habitats of humans, only cover 1% of the earth's surface, but they provide habitat for about 20% of all species worldwide [5,6]. Wetland ecosystem services are crucial for the survival of humans and a range of other animals, by supplying freshwater and aquatic products, flood control, climate regulation, biodiversity protection, etc. [7]. It has been estimated that lake wetland ecosystem services provide 23.2% of the total value of global ecosystem services [8].

Numerous scholars have assessed wetland ecosystems at different geographical locations and scales. For instance, Great Lakes Restoration Initiative Task Force systematically assessed the ecosystem services of the Great Lakes of North America in 2010, and showed that a number of management and engineering measures improved these services [9]. Wong et al. (2016) tested a new approach to evaluate water storage and local climate regulation of a green infrastructure project on the Yong ding River in Beijing, China, showing that the lakes and wetlands have increased local evapotranspiration on the Yongding Corridor by approximately 0.71 mm per day [10]. In another study, Bikash et al. (2015) pointed out that the economic values of the selective ecosystem services of the Koshi Tappu Wildlife Reserve, Nepal, equaled \$16 million per year [11]. Around the world, several countries, including Canada and China, have evaluated their wetland ecosystem services, but no universal evaluation index system of lake wetland ecosystem services has been established, and there has also been no evaluation of the effects of wetland ecosystem services on human health [12–17].

There is a lack of regional validation of parameters in the use of evaluation models, with only a few evaluations on cultural or support services in wetland ecosystem services, which makes the selection of indicators one-sided. In this sense, it is necessary to further unify the methods, indicators, and other aspects of ecosystem services of different scales and landscape types. At present, numerous studies are available regarding the quantitative evaluation, model simulation, and scenario analysis of ecosystem services, while using different types of ArcGIS (Aeronautical Reconnaissance Coverage Geographic Information System)-based support tools to fully integrate geography, ecology, economics, etc. The number and types of models for ecosystem services assessment have dramatically increased over the past decade. IMAGE (Integrated Assessment of Global Environmental Change), SAORES (Spatial Assessment and Optimization Tool for Regional Ecosystem Services), ARIES (Artificial Intelligence for Ecosystem Services), Costing Nature, and LUCI (Land Utilization and Capability Indicator) are some examples. For instance, ARIES is a model built via an expert knowledge method, which might result in the one-sidedness of evaluation results due to experts' knowledge reserve and subjective factors. The models IMAGE, Costing Nature, and LUCI are based on a global data set for ecosystem service evaluation [18,19]. The SAORES model evaluates ecosystem services in the context of ecological restoration and the management of the Loess Plateau in China, ignoring the accuracy verification and uncertainty analysis of service evaluation [20]. The InVEST (Integrated Valuation of Ecosystem Services and Trade-offs) model inputs the relevant data of the study area and reflects the impacts of ecosystem structure and function on service flow and value through the simplification of biophysical processes [21]. The multi-scale data input and the result output have made the InVEST model a comprehensive ecosystem service evaluation and balance model that is widely used in 20 countries and regions, for example in the assessment of freshwater ecosystem services in the Tualatin and Yamhill basins, the Chubute River valley in southern Argentina, the Beijing mountain forest, etc. [22–25].

With the rapid development of social economy and the expansion of global population, the demand for land is constantly increasing in order to survive and obtain economic benefits. Consequently, wetland resources are under unprecedented pressure and are decreasing at the global level [26,27]. The Dongting Lake Wetland is the largest lake that leads to the Yangtze River and the second largest freshwater lake in China. It provides favorable conditions for the breeding of snails and it is a severe epidemic area of *schistosomiasis*, seriously endangering the health of the local people. *Schistosomiasis* is a parasitic disease that affects both humans and animals and, globally, it is the second largest tropical parasitic disease after malaria [28,29]. It is endemic in about 78 countries in the tropical and subtropical regions of the world [30]. More than 95% of the *Schistosomiasis* in China are distributed in the Yangtze River basin tidal land, which seriously endangers wetland ecosystems. *Oncomelania hupensis*, which is a tropical freshwater snail, is the only intermediate host of *Schistosoma japonicum* and it can live in moist, shaded, and alternate land and water wetland environments. Thus, snail eradication is an effective measure for blocking the transmission of *schistosomiasis* [31,32]. At present, there are various technical

systems for snail control, including drug snail control, agricultural engineering measures, and water conservancy projects. The impermanence of snail control effectiveness, the higher costs on capital, and the labor-intensity are substantial disadvantages between water conservancy projects and ecological environmental protection. The 65th World Health Assembly has passed Resolution WHA65.21, which proposes eliminating *schistosomiasis*, a neglected tropical disease, in low-transmission areas of the world. Previous research has shown that the *Schistosomiasis* Prevented Forestry Project is an important measure for the comprehensive control of *schistosomiasis* in China, with the advantages of long-term effects and prevention [33–37]. The project is based on the construction of forests for snail control and *schistosomiasis* prevention, which aims to change the environment for snail breeding, inhibit the growth and development of snails, isolate the source of infection, and control the epidemic of *schistosomiasis* [38].

Recently, China has built a total of 5189 km2 of snail control and *schistosomiasis* prevention forests, decreasing the snail density and infection rate of human by 89.8% and 51.0%, respectively [39]. The construction of these forests, which also provide a variety of other ecosystem services, such as carbon fixation, water and soil conservation, and biodiversity protection, has reached 1640 km2, in 2015, along the Dongting Lake Wetland beach. However, such ecosystem services, especially snail control and *schistosomiasis* prevention, have been paid less attention in previous studies on the Dongting Lake Wetland area (e.g., Mao et al., 2007, where some ecosystem services, such as agricultural and fishery productions, tourism and leisure, climate regulation, and flood storage were investigated [40]). On the other hand, most studies regarding the snail control and *schistosomiasis* prevention in the literature have mainly focused on this issue in terms of mechanism and construction technology [39,41]. This study focuses on the ecosystem services in the Dongting Lake area that resulted from the construction of disease prevention forests, including snail control and *schistosomiasis* prevention, water yield, soil conservation, and carbon storage; and, evaluates them in a period of 10 years between 2005 and 2015. For this purpose, ArcGIS 10.2 and InVEST models (with the aid of market value and shadow engineering methods for service value estimations) are applied to (1) study the temporal and spatial variations of these ecosystem services in the study area; (2) explore their variations between different land use types; and, (3) estimate their value in 2005, 2010, and 2015.

#### **2. Study Area and Data**

#### *2.1. Study Area*

Dongting Lake is the second largest freshwater lake in China, which covers an area of 18,780 km2 and combines Xiangjiang, Zijiang, Yuanjiang, Lijiang, and other water-flooded lakes in the Yangtze River floods during the flood season. It is located between N 27◦39 ~30◦25 and E 111◦19 ~113◦34 , with an average elevation of 33.5 m above sea level. Average annual temperature in the lake is 17.0 ◦C, with an average annual precipitation of 1100–1400 mm. The soil pH ranges between 7 and 8.4. From June to August, the area is generally flooded; the difference between the highest water level in the wet season and the lowest one in the dry season in winter can reach 8–10 m, which makes the Dongting Lake Wetland a *schistosomiasis* epidemic area [42]. Figure 1 shows the geographic position of this area, which includes 17 counties and cities, such as Yueyang City, and Xiangyin, Huarong, and Li Counties in Hunan Province, covering an area of 3.17 <sup>×</sup> <sup>10</sup><sup>6</sup> km2, accounting for 16.5% of the area of the province.

By the end of 2010, the area of snail control and *schistosomiasis* prevention forests totaled 701.86 km2, while, in 2015, this area amounted to 1334.66 km2, which greatly inhibited snail numbers and decreased the number of *schistosomiasis* patients. The proportion of arable land in the Dongting Lake Wetland accounts for more than 45%, followed by woodland and water area. The proportions of residential, unutilized, and grassland areas are comparatively small and below 5%. From 2005 to 2015, arable land, woodland, grassland, and water area decreased at different degrees, i.e., by 1% (313.90 km2), 0.35% (111.09 km2), 0.15% (45.99 km2), and 0.72% (228.25 km2), respectively. In contrast, the proportions of

residential and unutilized lands continued to increase by 1.34% (424.26 km2) and 0.88% (279.60 km2), respectively. Table 1 presents the land cover types and areas from 2005 to 2015.

**Figure 1.** Geographic location of the study region (**a**) in China; (**b**) in the vicinity of Hubei, Hunan and Jiangxi provinces; and (**c**) elevation map of the study area including Dongting lake.


**Table 1.** Land cover types of the Dongting Lake Wetland between 2005 and 2015 (×10<sup>2</sup> km2).

#### *2.2. Data Sources and Processing*

The land use maps (at a spatial resolution of 30 × 30 m) of the Dongting Lake Wetland from 2000 to 2015 were downloaded from the Resource and Environment Data Cloud platform [43]. Precipitation data of the Wetland during this period were acquired from the National Meteorological Information Center [44]. ArcGIS 10.2 conducted Kriging interpolation in the study area and the surrounding areas on precipitation data from 59 meteorological and hydrological stations and, subsequently, grid maps with a resolution of 30 m were generated for the related years. The Digital Elevation Model (DEM) data were generated in ArcGIS 10.2 by means of interpolation, cutting, and filling digital elevation data of Global Digital Elevation Model Version 2 (GDEM V2), with a high resolution of 30 m, which was sourced from Geospatial Data Cloud [45]. Soil data, including soil texture, topsoil organic carbon, and fractions of topsoil sand, silt, and clay were obtained from the Soil Map of China based on the Harmonized World Soil Database (version 1.1, with a 1:1 million map scale), which was created by the Cold and Arid Regions Sciences Data Center in 2009 [46], and ArcGIS 10.2 was used to obtain the soil grid maps with a resolution of 30 m. The carbon density values for each land use were collected from Li et al. (2003) and Wang et al. (2010) [47,48] based on the principle of similarity within the natural environment of the Dongting Lake Wetland.

#### **3. Methods**

The data were firstly processed by using the ArcGIS 10.2 to develop an input for the analysis of the ecosystem services in this section, as discussed in the previous section. The InVEST model is then used, as described in Sections 3.1–3.3, to estimate water yield (the difference between precipitation and actual evapotranspiration), soil conservation, and carbon storage. Moreover, the Raster Calculator of ArcGIS 10.2 is applied to evaluate the snail control and *schistosomiasis* prevention, as presented in Section 3.4. Besides, the market value and shadow engineering methods were employed to calculate the values of the services.

#### *3.1. Water Yield Model*

The InVEST water yield module is an estimation method that is based on the water balance method, where the water yield of a grid unit refers to precipitation minus actual evapotranspiration, and the generated water yield data are the total amount and depth of water yield [49]. The module formula is as follows:

$$\mathcal{Y}\_{\mathbf{x}\mathbf{j}} = \left(1 - \frac{AET\_{\mathbf{x}\mathbf{j}}}{P\_{\mathbf{x}}}\right) \cdot P\_{\mathbf{x}} \tag{1}$$

where *Yxj* and *AETxj* are the annual water yield and the actual evapotranspiration of the land cover type *j* in Grid *x* in mm; and, *Px* is the precipitation of Grid *x* in mm. The annual potential evaporation for the 87 meteorological stations in the Dongting Lake Wetland was estimated by the Penman–Monteith equation [50]. Subsequently, Kriging spatial interpolation was used to obtain a map of the annual average potential evapotranspiration of the Wetland. Besides, the market value method was employed to calculate the value of the water yield service.

#### *3.2. Soil Conservation*

The soil retention module, which was based on the Revised Universal Soil Loss Equation (RUSLE), was used to estimate the amount of soil erosion within the region. The formulas of the module are as follows:

$$A = R \times K \times LS \times \mathbb{C} \times P \tag{2}$$

$$RKLS = R \times K \times LS \tag{3}$$

$$\text{LSD} = \text{RKLS} - A \tag{4}$$

where *RKLS* is the total potential soil loss (t) per pixel in the original land cover without the application of factors *C* or *P* from the RUSLE (i.e., equivalent to the soil loss for bare soil); *A* is the estimated average soil loss (t) per pixel in the original land cover (t·km−2·yr<sup>−</sup>1); *SD* is the amount of soil retention (t); and, *<sup>R</sup>* is the rainfall erosivity (MJ·mm/(km2·h·yr)), which is calculated while using the formula of Wischmeier according to the average monthly and annual precipitations in the Dongting Lake Wetland [51]; *K* is the soil erodibility (t·km2·h/(km2·MJ·mm)) calculated while using the formula of Williams based on relevant soil data [52]; *LS* is the slope length gradient factor; *C* is the cover management factor obtained by the FAO [53]; and, *P* is the support practice factor that is available from the USDA Agriculture Handbook [54]. The value of soil conservation was estimated while using shadow engineering method.

#### *3.3. Carbon Storage*

The carbon module uses land use maps and the carbon density (i.e., above ground biomass, belowground biomass, soil, and dead biomass) of each land use type to estimate the amount of carbon that is currently stored in a landscape, as follows [55]:

$$\mathcal{C}\_t = \mathcal{C}\_{abvve} + \mathcal{C}\_{belav} + \mathcal{C}\_{dend} + \mathcal{C}\_{soil} \tag{5}$$

where *Ct* is the total carbon storage (t); and, *Cabove*, *Cbelow*, *Csoil*, and *Cdead* are the carbon storage (t) in aboveground biomass, belowground biomass, soil, and dead biomass, respectively. The value of carbon storage was calculated while using the market value method.

#### *3.4. Snail Control and Schistosomiasis Prevention*

The purpose of constructing snail control and *schistosomiasis* prevention forests has been to reduce the density and area of snails and effectively prevent and control the spread of *schistosomiasis*. Therefore, in this study, the snail-control function is evaluated by the reduction of the snail area and the reduction in the number of *schistosomiasis* patients. The Raster Calculator of ArcGIS 10.2 is used to analyze the quality and value of snail control and disease prevention, while using the following equations:

$$
\Delta\_A = A\_{2010} - A\_{2015} \tag{6}
$$

$$
\Delta\_N = N\_{2010} - N\_{2015} \tag{7}
$$

where Δ*<sup>A</sup>* and Δ*<sup>N</sup>* denote the reduced areas that are colonized by the snail and the reduction number of *schistosomiasis* patients, respectively; *A*<sup>2010</sup> and *N*<sup>2010</sup> are the snail area (104 m2) and the number of schistosomiasis patients in 2010; and, *A*<sup>2015</sup> and *N*<sup>2015</sup> are those quantities in 2015.

The market value and the shadow engineering methods, and then were used to calculate the value of snail control and *schistosomiasis* prevention. The value of snail and disease prevention function (*C*) is computed as the following:

$$\mathcal{C} = \mathbb{C}\_{\text{Small}} \cdot \Delta\_A + \mathbb{C}\_{\text{Patient}} \cdot \Delta\_N \tag{8}$$

where *CSnail* is the cost of snail control (in CNY/m2) and *CPatient* is the annual treatment cost per patient. Table 2 summarizes the sources used for the collection of parameters.


#### **Table 2.** Parameter collection methods.

#### **4. Results and Analysis**

In this section, the results of the present study on the variations of different ecosystem services in the Dongting Lake Wetland and an analysis of their values are presented.

#### *4.1. Temporal and Spatial Variations of the Ecosystem Services*

Figure 2 shows the results of the study on the temporal and spatial variations of the water yield, soil conservation, and carbon storage; and, Figures 3 and 4 present those variations in the reduced snail area and the number of schistosomiasis patients, in the Dongting Lake Wetland, from 2005 to 2015.

**Figure 2.** Distribution of changes in the ecosystem services in the Dongting Lake Wetland from 2005 to 2015.

**Figure 3.** Temporal and spatial variations of the reduced snail area from (**a**) 2005 to 2010; (**b**) 2010 to 2015; and (**c**) 2005 to 2015.

**Figure 4.** Temporal and spatial variations of the reduction in the number of *schistosomiasis* patients from (**a**) 2005 to 2010; (**b**) 2010 to 2015; and (**c**) 2005 to 2015.

#### 4.1.1. Variations of Water Yield

Land use and climatic changes both influence water yield. In 2005, 2010, and 2015, water yield in the Dongting Lake Wetland was 236.91 <sup>×</sup> 10<sup>8</sup> m3, 339.61 <sup>×</sup> 10<sup>8</sup> m3, and 302.80 <sup>×</sup> 108 m3, respectively. This shows that it increased from 2005 to 2010 by 43.3% and then decreased 10.8% between 2010 and 2015. The reason for these variations was that the amount of rainfall was higher than evapotranspiration

between 2005 and 2010, while the latter became larger from 2010 to 2015. However, the net change from 2005 to 2015 was an increase of 65.89 <sup>×</sup> 108 m3 (equivalent to 27.81%). Spatially, water yield during this period decreased in the southwest (Taoyuan and Hanshou Counties), with the largest reduction of 1313.0 m3/km2, while it increased in the northeast (Yueyang County, and Linxiang and Yueyang Cities), where the largest rise was 1586.8 m3/km2.

#### 4.1.2. Variations of Soil Conservation

In 2005, soil conservation was equivalent to 1257.23 <sup>×</sup> <sup>10</sup><sup>4</sup> t, while, in 2010 and 2015, it accounted for 1931.06 <sup>×</sup> 104 t and 1257.44 <sup>×</sup> 10<sup>4</sup> t, respectively, i.e., a total increase of 0.21 <sup>×</sup> 104 t equivalent to a growth of 0.02% from 2005 to 2015. However, the change was not uniform across the area, as in the northwest (Li County) and southwest (Taoyuan County) a decrease was observed with the largest value of 4569.7 t; while, in the northeast (Linxiang City, Yueyang County, and Yueyang City), the amount of conserved soil became larger (with a maximum change of 61,943 t). See Figure 2 for details.

A variety of factors, such as climate, topography, land use, and human activities, influence soil conservation. Among these, precipitation and rainfall erosivity, with temporal variations that were similar to the soil conservation (i.e., increase from 2005 to 2010 and then decrease from 2010 to 2015), were the most significant causes of the large variations in the amount of the conserved soil between 2005 and 2015.

#### 4.1.3. Variations of Carbon Storage

In 2005, 2010, and 2015, the total carbon storage was 3.57 <sup>×</sup> 108 t, 3.63 <sup>×</sup> 10<sup>8</sup> t, and 3.56 <sup>×</sup> 108 t, respectively. This shows a total decrease of 0.01 <sup>×</sup> <sup>10</sup><sup>8</sup> t (i.e., 0.28%) from 2005 to 2015. Spatially, during this period, the carbon storage mainly changed in the northwest (Li County), northeast (Yueyang County), and southeast (Wangcheng County), with a maximum reduction of 26.2215 t; while, in the areas in the middle (Huarong county and Yuanjiang City), it increased with a maximum change of 26.2215 t (Figure 2). Temporal and spatial changes in the carbon storage were both associated with the changes in the land use. For example, higher variations in the carbon storage were observed in the regions with higher changes of land use.

#### 4.1.4. Variation of Snail Control and Schistosomiasis Prevention

The construction of snail control and *schistosomiasis* prevention forests in the Dongting Lake Wetland has considerably changed the environmental conditions for snail breeding. Figure 3 shows the temporal and spatial reductions of snail areas in the study area between 2005 and 2015. In 2005, snails colonized 176,042.5 <sup>×</sup> 104 m<sup>2</sup> of the area. Subsequently, from 2005 to 2015, 10 counties experienced a reduction, with the total amount of 326.01 <sup>×</sup> 104 m<sup>2</sup> (i.e., <sup>−</sup>0.19%), while an increase of snail area was observed in seven other counties. These variations were as follows. In Taoyuan, Anxiang, Wangcheng, Huarong, Hanshou, and Xiangyin Counties, as well as Yueyang, Yuanjiang, Jin, and Linxiang Cities, there was a reduction in the snail area with a maximum of 820.37 <sup>×</sup> 104 <sup>m</sup><sup>2</sup> in Taoyuan County; while, in Nan, Yueyang, and Li Counties and Yiyang, Miluo, Linli, and Changde Cities, the area that was colonized by snails became greater with a maximum of 700.5 <sup>×</sup> 104 m2 in Li City. The increase in the snail area in the latter, in spite of the presence of the constructed snail control and *schistosomiasis* prevention forests, was that these areas are located along the *Oncomelania hupensis* snail-ridden regions outside embankments of Yangtze river and Dongting Lake, where the infectious sources are various, large in number, and difficult to manage, thus a high risk of schistosomiasis still exists [57].

Except in Linli County, where the number of *schistosomiasis* patients increased from 2005 to 2015 (386 people), this number decreased in all other counties and cities of the study area. In 2005, the total number of *schistosomiasis* patients in the Dongting Lake Wetland area was 210,651, which then reduced by 182,829 (86.8%) from 2005 to 2015. The largest reduction was observed in Yuanjiang City with 27,586 patients, followed by Nan County, while the smallest decrease was found in Wangcheng County with 154 patients. The increase in the number of patients in Linli County was mainly due to

the fact that, in this county, there is a large number of cattle and sheep, which increases the human and animal infection rate [58]. Figure 4 shows these temporal and spatial variations in the study area between 2005 and 2015.

#### *4.2. Variations in Ecosystem Services between Di*ff*erent Land Use Types*

Different types of land use showed different changes in the water yield. The amount of water yield per km<sup>2</sup> in unutilized, residential, arable, forest, and grassland areas had a growth of 63.83%, 32.86%, 28.19%, 23.25%, and 25.79%, respectively (Figure 5a). Moreover, it can be seen from this figure that the order of water yield per km2 in different land use types was as follows: woodland > arable land > grassland > residential land > unutilized land > water area.

**Figure 5.** Ecosystem services per km<sup>2</sup> of different land cover types in the Dongting Lake Wetland from 2005 to 2015: the amounts of (**a**) water yield; (**b**) soil conservation; and (**c**) carbon storage.

The change in the total amount of conserved soil was different every five years. As shown in Section 4.1.2, it increased by 53.60% in the first five years of the study period, while it decreased by 34.88% in the second 5 years. Soil conservation per km2 of different land use types from 2005 to 2015 is shown in Figure 5b. Among them, the residential, unutilized and grassland areas experienced an increase of soil conservation per km<sup>2</sup> by 38.74%, 37.40%, and 12.76%, respectively, while the others had a smaller increase. In summary, the soil conservation per km2 in different land use types between 2005 and 2015 followed this order: grassland > woodland > water area > arable land > residential land > unutilized land.

As shown in Section 4.1.3, from 2005 to 2010, the total amount of carbon storage increased by 1.68%, and then a decrease of 1.93% was observed from 2010 to 2015. Figure 5c presents the changes in the carbon storage capacity per km<sup>2</sup> of different land use types for every five years during the study period. As can be seen, there has been an increase in unutilized land and grassland by 3.94% and 1.14%, respectively; and, a decrease of 4.87%, 0.15%, and 0.03%, respectively, in the residential, arable, and woodland areas. According to the study of Chen et al., (2012) [59] on the carbon storage in Dongting Lake area from 1979 to 2004, soil organic carbon was largely lost when grassland was converted into woodland due to the establishment of forests in the Wetland area. According to the present results, the

order of the carbon storage was as follows: grassland > woodland > unutilized land > residential land > arable land.

#### *4.3. Analysis of the Value of Ecosystem Services*

Figure 6 presents the results of the analysis that was carried out on the value of the ecosystem services in the Dongting Lake Wetland while using the market value and shadow engineering methods. The value of these services was estimated to be 1249.84 <sup>×</sup> 108 CNY in 2005, of which water yield, carbon storage, soil conservation, and snail control and *schistosomiasis* prevention accounted for 85.49%, 13.87%, 0.64%, and 0.00%, respectively. In 2010, this value increased to 1723.85 <sup>×</sup> 108 CNY, with the rise in the value of water yield, soil conservation and snail control, and *schistosomiasis* prevention by 3.36%, 0.07%, and 0.23%, respectively, and a reduction of 3.67% in carbon storage. In 2015, the overall ecosystem value reduced to 1552 <sup>×</sup> 108 CNY. This was associated with decreases of 0.88% and 0.20% in the water yield and soil conservation, respectively, and growth of 0.93% and 0.15% in the values of carbon storage and snail control and *schistosomiasis* prevention.

**Figure 6.** Values of different ecosystem services of the Dongting Lake Wetland in (**a**) 2005; (**b**) 2010; and (**c**) 2015.

Besides, it was observed that the values of different land use types per km2 followed the order of snail control and *schistosomiasis* prevention forest > forest (non-snail control and *schistosomiasis* prevention) > arable land > residential land > unutilized land > water area. This can be seen in Figure 7, where the values per km2 of different land cover types are presented for years 2005, 2010, and 2015.

**Figure 7.** Values per km<sup>2</sup> of different land cover types.

In 2010, the value of the ecosystem service of snail control and *schistosomiasis* prevention forests in the Dongting Lake Wetland area was 54.33 <sup>×</sup> 108 CNY, which corresponded to the water yield, soil conservation, carbon storage, and snail control, and *schistosomiasis* prevention by 74.25%, 1.10%, 17.31%, and 7.34%, respectively. In 2015, this value reached 89.11 <sup>×</sup> 108 CNY. Only the carbon storage value increased from 2010 to 2015 (by 2.33%), but the value of the other factors reduced, as seen in Figure 8.

**Figure 8.** Values of different ecosystem services of snail control and *schistosomiasis* prevention forests in the Dongting Lake Wetland area in (**a**) 2010; and (**b**) 2015.

#### **5. Discussion**

#### *5.1. Factors A*ff*ecting the Ecosystem Services*

Previous studies have found that land use change, precipitation, evapotranspiration, and rainfall erosivity are the main mechanisms that affect ecosystem services [60]. We have carried out a correlation analysis of the ecosystem services with evapotranspiration, precipitation, soil erodibility, and rainfall erosivity while using Statistical Product and Service Solutions (SPSS) to investigate the effects of such factors in the present study. Table 3 presents the result, where it can be seen that, from 2005 to 2015, the ecosystem services show different degrees of correlation with these factors. Specifically, all of the factors significantly positively influenced water yield, soil conservation, and carbon storage; and, the rainfall erosivity positively influenced the reduced snail area at a moderate degree.


**Table 3.** Correlation analysis of ecosystem services and various influencing factors during 2005–2015.

Asterisks indicate significant differences at \*\* *p* < 0.01.

#### *5.2. Selection of Appropriate Evaluation Indexes and Methods*

Wetlands provide a variety of ecosystem services. The evaluation indexes of wetland ecosystem services are different due to different characteristics of wetlands. For example, Schallenberg et al. (2013) conducted ecosystem services evaluation on eight lakes in New Zealand, with an area of over 100 km2 including water supply, fish supply, water purification, flood control, lakeside protection, climate regulation, and recreation, and discovered that they were in the trend of degradation [61]. Zhang et al. (2014) calculated the ecological service value of lakes and wetlands in China from three aspects of biodiversity, water quality, and economic indices. The results showed that the values of lake and marsh wetlands ecosystem services accounted for 54.64% of the total value of China's natural grassland ecosystem services and 30.34% of the total value of forest ecosystem services, respectively [62]. Sun et al. (2017) compared the different results of food security and biodiversity between China's Poyang Lake and Bangladesh's Tanguar Haor wetlands and showed that the trends in food security and biodiversity are declining [63]. Wondie (2018) selected nine wetlands in the region of the Lake Tana (Ethiopia) and estimated the value of wetland ecosystem services in terms of supply, regulation, culture, and support services. They showed that the wetland utilization planning, investment intervention, public land ownership, and other policy decision-making supports are poor, thus resulting in a major challenge to the sustainable use of wetland resources [64]. Mao et al. (2007) estimated the value

of the ecological services in the Dongting Lake Wetland to be 411.10 <sup>×</sup> 108 CNY using the market value, shadow engineering, carbon tax, and alternative expenses methods, but without considering the effect of snail control and *schistosomiasis* prevention forests [40]. We observed that in 2015, this value was 1552 <sup>×</sup> 108 CNY, and the snail control and *schistosomiasis* prevention forests had the highest ecosystem service value among the different land use types. This is mainly due to the fact that such forests provide various benefits, such as water yield, soil conservation, and carbon storage. However, more importantly, they control and prevent *schistosomiasis*, which thereby reduces the cost of snail eradication and *schistosomiasis* treatment and decreases the income loss of *schistosomiasis* patients.

In summary, different assessment indicators are adopted for wetlands in different areas, and quite different results are often obtained. It can be concluded that there is no consensus on a singular best approach for the assessment of wetland ecosystem services, but the selection of appropriate assessment indicators and methods for a certain wetland is heavily dependent on the characteristics of the area where the wetland is located.

#### *5.3. Limitation of the Present Study*

Our study has some limitations, as follows. (1) This study only evaluated water yield, soil conservation, carbon storage, and snail control and *schistosomiasis* prevention functions in the evaluation of the ecosystem services of snail control and *schistosomiasis* prevention forests. However, other ecosystem services, such as the provision of wood and forest products, climate regulation, and air purification, which may have significant effects, were neglected. (2) The dynamic impact of afforestation years on the ecosystem services, such as the difference between an old-growth forest and a newly planted forest, was disregarded.

#### **6. Conclusions**

In this study, water yield, soil conservation, carbon storage, and snail control and *schistosomiasis* prevention were evaluated in the Dongting Lake Wetland area from 2005 to 2015. Their temporal and spatial changes were studied in detail, where different trends and variations were observed. Different land use types also showed various changes in these services. Subsequently, the value of the ecosystem services was estimated, and a total increase of 24.18% was observed over the 10-year period, where the water yield showed the highest value, which accounted for more than 70% of the total value of the four evaluated ecosystem services; while, the value of the snail control and *schistosomiasis* prevention was about 5% of it. Moreover, through a correlation analysis, it was shown that evapotranspiration, precipitation, soil erodibility, and rainfall erosivity significantly positively influenced some of the ecosystem services in the study area.

There is no singular best approach for the assessment of ecosystem services of wetlands, but the appropriate methods should be selected based on the characteristics of the area where the wetland is located. The potential improvements in the selection of evaluation index system and methods for a similar case to our study include considering the dynamic impact of afforestation years on the ecosystem services and taking the effect of other ecosystem services, such as the provision of wood and forest products, climate regulation, and air purification into account.

**Author Contributions:** Conceptualization: L.M., R.S., D.P.; methodology, software, validation and writing original draft preparation: L.M.; writing—review, editing and revision: L.M., E.K., Y.Z., D.P., J.Z., K.Z.; funding acquisition: J.Z.; project administration: J.Z.; supervision: Q.S.

**Funding:** The National Science and Technology Support Program of China (grant No. 2015BAD07B07), and the National Natural Science Foundation of China (NSFC) (grant No. 41071334) supported this research.

**Acknowledgments:** We would like to thank Yilin Zhuang for her support and suggestions in the preparation of the manuscript.

**Conflicts of Interest:** Authors declare there is 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* **Study on Multi-Scale Coupled Ecological Dispatching Model Based on the Decomposition-Coordination Principle**

#### **Tao Zhou \*, Zengchuan Dong, Wenzhuo Wang, Rensheng Shi, Xiaoqi Gao and Zhihong Huang**

Hydrology and Water Resources College, Hohai University, Nanjing 210098, China **\*** Correspondence: 160401010016@hhu.edu.cn

Received: 9 May 2019; Accepted: 9 July 2019; Published: 12 July 2019

**Abstract:** Studies on environmental flow have developed into a flow management strategy that includes flow magnitude, duration, frequency, and timing from a flat line minimum flow requirement. Furthermore, it has been suggested that the degree of hydrologic alteration be employed as an evaluation method of river ecological health. However, few studies have used it as an objective function of the deterministic reservoir optimal dispatching model. In this work, a multi-scale coupled ecological dispatching model was built, based on the decomposition-coordination principle, and considers multi-scale features of ecological water demand. It is composed of both small-scale model and large-scale model components. The small-scale model uses a daily scale and is formulated to minimize the degree of hydrologic alteration. The large-scale model uses a monthly scale and is formulated to minimize the uneven distribution of water resources. In order to avoid dimensionality, the decomposition coordination algorithm is utilized for the coordination among subsystems; and the adaptive genetic algorithm (AGA) is utilized for the solution of subsystems. The entire model—which is in effect a large, complex system—was divided into several subsystems by time and space. The subsystems, which include large-scale and small-scale subsystems, were correlated by coordinating variables. The lower reaches of the Yellow River were selected as the study area. The calculation results show that the degree of hydrologic alteration of small-scale ecological flow regimes and the daily stream flow can be obtained by the model. Furthermore, the model demonstrates the impact of considering the degree of hydrologic alteration on the reliability of water supply. Thus, we conclude that the operation rules extracted from the calculation results of the model contain more serviceable information than that provided by other models thus far. However, model optimization results were compared with results from the POF approach and current scheduling. The comparison shows that further reduction in hydrologic alteration is possible and there are still inherent limitations within the model that need to be resolved.

**Keywords:** ecological operation; multi-scale; decomposition-coordination; hydrologic alterations

#### **1. Introduction**

Reservoir ecological operation has become increasingly more significant and has recently drawn a lot of attention from research scientists in numerous scientific disciplines. By 2005, there were >45,000 reservoir dams higher than 15 m throughout the world, which are collectively capable of storing ~15% of annual global runoff [1]. While these reservoirs make a sizeable contribution to developing the economy, they are also the source of numerous ecological problems. The ultimate goal of ecological operation is to integrate ecological factors into standard operation objectives, such that ecological protection is given the same consideration as increasing profit margins. Essentially, the target of ecological operation is to mediate the conflict between society's demand and the river ecosystem's needs [2]. As early as 1971, Schlueter put forward that reservoir management should maintain river diversity [3]. Following that pronouncement, a number of river ecological water demand calculation methods were developed—i.e., the Tennant Method, Wetted Perimeter Method, and IFIM Method, etc. [4]. Since then, increasing discharge from reservoirs to reserve enough ecological flow has become a common method of ecological operation [5–7].

Progressively more in-depth research on the riverine ecosystems has widened our understanding of the river system itself and the ecosystem it supplies. Thus, river ecological water demand no longer exclusively concerns the reserved water downstream. Many scholars have reached a broad consensus that riverine ecosystems adapt to the natural river flow pattern [8], and every river has its own characteristic flow pattern and corresponding biological community [9]. Thus, any flow pattern change—such as flood runoff, flow velocity, arrival time, and flood pulse frequency, etc., can cause a series of ecological consequences. As a result of these insights, ecological water requirements have been transformed from merely considering reserved flow [10,11] to considering the full extent of natural flow to which riverine biological communities must adapt [11–13]. Based on the above consensus, Richter et al., 1997 suggested implementing the Range of Variability Approach (RVA), which uses the Indicator of Hydrologic Alteration (IHA) to quantify how much the reservoir operation alters the riverine ecosystem. Homa et al. [11] presented the concept of ecological deficit and proposed minimizing ecological deficit as an objective of reservoir operation. Richter [14] offered a definition of sustainable water management called "Sustainability Boundary Approach (SBA)". Francisco et al. used SBA to assess and manage three major river basins in Spain with good results [15,16]. Typical reservoir operating rule curves (RORCs) [17] were improved by using RVA as evaluation indicators [18,19].

While it is now recognized that a comprehensive evaluation of ecological disruption requires an investigation over the entire river system, the water need in river ecosystems must also be considered on multiple scales. The minimum ecological discharge is mainly reflected on a monthly scale, but pulse flow and ecological flow velocity are reflected on smaller scales. As such, studies that are restricted to any single scale are inherently missing important ecological information. That being said, the large scales used for current research on reservoir ecological operation are too high-level to describe ecosystem demand with more than an overview, as flow pattern and flood pulse, etc., are too specific to be considered.

Given that the ecological water demand is present on multiple scales, any ecological dispatching model needs to include various scales. Steinschneider et al. [20] developed a large-scale optimization model to evaluate the contribution of reservoir management measures to ecological benefits; and Tsai et al. [21] proposed a new hybrid method to quantify river ecosystem demand based on artificial intelligence. Clearly, the importance of a large-scale optimization model and small-scale characteristic flow regimes has been recognized, but few studies couple them together in a feasible and practical way. Although the multi-scale dispatching model covers more ecological information, the number of model decision variables dramatically increases accordingly. For example, the number of decision variables in this study is nearly 3000. Increasing the number of model decision variables and the non-linear forms of hydrologic alteration makes the overall system model difficult to solve and even harder to understand. Furthermore, current studies usually employ larger step sizes (such as ten-day period) and take scheduling rules parameters as decision variables [17], but they fail to acquire all the characteristics of daily flows. Finally, the large-scale dispatching model is unable to be learned by machine-learning algorithms. In order to overcome these challenges, a hierarchical control method of the decomposition-coordination principle was applied in this study. The large system decomposition coordination theory is a new branch of science that has developed rapidly since the 1970s. In its most fundamental form, the theory proposes decomposing a complicated system into several independent subsystems and handling the correlation between subsystems by coordinators. In 1984, Adigüzel [22] decomposed a complicated system into several subsystems according to the objective function based on the decomposition-coordination principle. Currently, the decomposition-coordination principle is one of the most commonly used methods for solving complex systems. Over time, as the decomposition-coordination principle was more frequently applied, a

multi-layer hierarchical control system was established. In the multi-layer hierarchical control system, the superior subsystem result is not only the boundary but is also harmonized with the result of the subordinate subsystem.

As stated above, researchers acknowledge the multi-scale features of ecological water demand [23,24]; yet few studies have been performed on the multi-scale coupled ecological dispatching model. Thus, in order to obtain a more reasonable deterministic optimization model, a multi-scale coupled ecological dispatching model based on the decomposition-coordination principle was built in this study. The scheduling model consists of large-scale and small-scale models, which use the month and day, respectively, as the calculation scale. The large-scale model was implemented for solving problems associated with water supply and river ecological base flow; while the small-scale model was used for solving problems related to extremely low flows, flood pulses, and ecological hydraulic parameters. Due to the introduction of RVA, disturbances are kept within reasonable bounds, which is beneficial to the river ecosystems. Furthermore, the simulative daily flows can be applied to the training of ecological regulations by a machine learning algorithm, which is important to the implicit stochastic optimization. Investigations that take scheduling rules parameters as decision variables or use larger step sizes cannot obtain valid daily flow data. The multi-scale coupled model also provides a framework for the formulation of multi-scale scheduling rules. This model can help decision-makers to determine whether there is still room for improvement in ecological scheduling and look for ways to improve.

#### **2. Materials and Methods**

#### *2.1. Generation and Evaluation of Characteristic Flow Regimes*

The importance of natural flow regime to the riverine ecosystems has been widely accepted. Richter et al. [4,25] put forward IHA and RVA indicators to quantify the degree of flow regime alteration. Consequently, the effect of reservoir ecological operation can be quantified. Mathews and Richter [26] suggested that flows can be divided into five categories: small floods (flows ≥ the overbank flow, but <10 year flood), high-flow pulses (flows < the overbank flow, but > seasonal base flow), low flows (base flows in different months), extreme low flows (flows ≤ 95th percentile flow) and large floods (flows ≥ 10 year flood). Yin et al. [19] combined small floods and large floods into one flood category. Furthermore, to consider additional physical mechanisms, we take into account the hydraulic parameters that do not belong to the IHA in this study.

In order to build a comprehensive multi-scale coupled ecological dispatching model, the flow regime alteration level must be quantified, and the approaches for generating the characteristic flow regime must also be explained. In this study, all environmental flow components were accounted for by considering four kinds of characteristic flow regimes: base flows, extreme low flows, flood pulse with high flow pulses, and ecological hydraulic parameters. The generation mechanisms for these characteristic flow regimes are based on their ecological significance.

#### 2.1.1. Base Flows

Ecological functions of base flows include maintaining suitable water quality, sustaining hydraulic habitats for aquatic organisms, and supporting hyporheic organisms [19]. Thus, maintaining base flow in a river is very important for sustaining the riverine ecosystem. Based on the ecological functions of base flow, the full spectrum of base flow magnitudes was selected as evaluation indicators in the multi-scale coupled ecological dispatching model.

Depending on the dates in question, a method (e.g., Tennant or 7Q10 methods, etc.) was selected and implemented to calculate the base flows. The base flows were then used as constraint conditions in the large-scale scheduling model, which has an account step of 1 month. Merging the large- and small-scale factors guaranteed that monthly discharge ≥ base flows.

With respect to the model's objective functions, RVA was used to quantify the degree of deviation between the discharged flow's base flow and the natural flow regime. In this study, the daily flow from 1958 to 1977 were used as the reference series—as the values closely reflect natural runoff; 2007–2014 was selected as the scheduling period. In order to avoid ecological flow flattening, the full spectrum of base flows within the reference period was divided into three target ranges using the 25th and 75th percentiles as parameter values. The parameters in every target range were preserved. The corresponding target ranges were integrated into the objective functions and mathematically expressed as follows:

$$D\_i = \sum\_{La=1,2,3} \left| \frac{N\_i(La) - N\_\ell(La)}{N\_\ell(La)} \right| \tag{1}$$

where *La* = 1, 2, 3 represents the high, middle and low target ranges, respectively; *Di* refers to the degree of change of the *i*-th index; *Ni*(*La*) refers to the number of years in which the i-th index is observed in the corresponding La target ranges; *Ni*(*La*) refers to the number of years in which the i-th index is expected to be observed in the corresponding *La* target ranges which were calculated according to the reference series.

As the number of parameters from the simulation period that are observed in the target ranges increases, the degree of hydrologic alteration decreases in parallel.

#### 2.1.2. Extreme Low Flows

Extreme low flows concentrate prey species for predators, provide habitat for terrestrial animals, and purge exotic species. At the same time, the extremely low flows can reduce the connectivity of rivers and confine the movement of certain aquatic organisms.

For ecological reasons, the magnitude, frequency, duration, and timing of extremely low flows were taken into consideration in the small-scale scheduling model. As with base flows, the full spectrum of the parameters (the magnitude, frequency, duration, and timing) was grouped into three target ranges using the 25th and 75th percentiles as parameter values. The mathematical expression is identical to that shown in Formula (1). Just as with the base flows, the more parameters from the simulation period that are observed in target ranges, the lower the degree of hydrologic alteration.

#### 2.1.3. Flood Pulses and High Flow Pulses

Junk [27] introduced the Flood Pulse Concept (FPC). Flood pulse is the primary driving force behind biological survival, productivity, and interaction of the river-flood plain system. The flood pulse's activity has the following ecological significance: (1) It promotes energy exchange and material circulation between aquatic species and terrestrial species; (2) it improves the dynamic connectivity of river systems, which provide shelter and habitat for different species; and (3) it transmits abundant and intense information to species, which is important for spawning, hatching, growth, refuge, and other migration activities.

Flood pulse and high flow pulse generation mechanisms are similar; they differ in threshold values. Using flood pulse as an example, in order to represent its ecological purposes, the magnitude, duration, timing, and frequency are considered in the small-scale scheduling model. Thus, the flood pulse generation mechanism is as follows:

When the control section flow exceeds a certain limit, a flood pulse incident is triggered.

$$\mathbf{Q}\_{\rm T} > \mathbf{T} \tag{2}$$

where T is a pre-set limit value that depends on the stage-discharge relationship and the control section shape; and QT is the control section flow.

When the control section flow decreases below a certain limit, the flood pulse incident ends. In the small-scale scheduling model, a flood pulse is induced by artificial floods due to the release flow from reservoirs. Thus, reservoirs need to control the release to compensate for the flood pulse as follows.

$$\mathbf{Q}\_{\rm P} = \mathbf{Q}\_{\rm R} - \mathbf{Q}\_{\rm T} \tag{3}$$

where QP represents the released flow from reservoirs; QR represents the control section target flow; and QT represents the control section flow.

Since flood pulse incidents usually last for a few days or less, conventional reservoir scheduling models with 1 month account steps, cannot incorporate the requirement. Thus, the magnitude, frequency, timing, and duration of all flood pulse incidents during the ecological critical period were evaluated in the small-scale scheduling model. The mathematical expression is the same as that shown in Formula (1).

#### 2.1.4. Ecological Hydraulic Parameters

Using hydraulic methods to calculate ecological water demand is both prominent and important. Hydraulic methods are based on hydrodynamic parameters as well as specific biological criteria. Conventional hydraulic methods include the wet perimeter method and R2-CROSS [7,28–30]. Liu et al. [31] introduced the concept of ecological flow velocity, which combines hydrodynamic parameters and specific ecological targets.

Ecological flow velocity differs from the above-mentioned characteristic flow regimes in that the ecological hydraulic parameters reference standard is based on both the natural flow regime and species requirements. Because ecological flow velocity reflects the ecological target requirement, it is expected to vary within a certain range. In addition, specific water depths—which will differ depending on the local channel morphology and the local ecology—are essential for maintaining the basic ecological function of the ecosystem. In order to fully incorporate the biological requirements and hydrodynamic parameters, this study accounts for the appropriate hydraulic parameters, such as the depth and flow velocity of specific species in different growth cycles. The selected species primarily represent those that dominate the riverine ecosystem. For example, the growth cycle of *Gymnocypris eckloni herzenstein*—a keystone species found upstream of the Yellow River—was divided into three stages: spawn, brood, and growth. Liu [32] determined that maintaining a larger flow velocity during the spawning period stimulates spawning. In contrast, slow to moderate flow velocity is optimal for the brood or growth stage. In terms of depth, Xu [33] established that fish require a minimum water depth more than three times of their length to survive (See Table 1).

**Table 1.** Gymnocypris eckloni herzenstein.


Table 1 shows ecological parameters of the *Gymnocypris eckloni herzenstein*. This data is from the Advanced Aquatic Biology [32].

Flow was related to the hydraulic parameters via an open channel uniform flow theory. The conventional reservoir scheduling models were unable to meet the requirement because the hydraulic parameters are sensitive to flow rate. Therefore, the small-scale scheduling model uses a day scale as the account step. The ecological hydraulic parameters evaluation function is the same as that shown in Formula (1). However, in this case, the target range is based on the natural flow and species requirements. The more hydraulic parameters, from within the corresponding natural flow values that are observed in the suitable ranges, the lower the degree of hydrologic alteration that will be obtained.

The descriptive parameters of the four discussed characteristic flow regimes are summarized in Table 2:



Table 2 shows descriptive parameters of characteristic flow regimes and which model will handle them.

In this study, the characteristic flow regimes' degree of alteration was averaged and used as the degree of hydrologic alteration.

#### *2.2. Construction of the Multi-Scale Coupled Ecological Dispatching Model*

Two key problems with the multi-scale coupled model include: (1) increased decision variables and (2) coordination between the different scale models. In order to resolve these issues, the decomposition-coordination principle and multiple hierarchical control system, which is based on the theory of large-scale systems, were applied to this work.

#### 2.2.1. The Space and Time Decomposition-Coordination Method

To resolve the issue of increased decision variables due to the introduction of small-scale scheduling models, the whole basin scheduling model for the entire scheduling period was divided into numerous subsystems, each of which includes part basin and part scheduling period. Because the whole basin scheduling model is regarded as a large complex system, the connection between the overall system and subsystems is reflected by the subsystems' objective functions, which include Lagrange multipliers. In this way, the large system optimization problem is transformed into several optimization problems of the subsystem objective functions. This method greatly reduces decision variables while retaining the model's optimal conditions [34]. With respect to space, the large-scale system is usually divided by reservoir control basins. In terms of time, the entire scheduling period is divided into small periods according to the model requirements. In general, subsystems coordinate with each other via subsystem hydraulic connections and cooperative work between subsystems—i.e., the interconnection constraints. The subsystem is mathematically expressed as follows:

$$L = \sum\_{i=1}^{n} L\_i = \sum\_{i=1}^{n} f\_i(\mathbf{x}\_i) + \lambda \left[ \sum\_{i=1}^{n} g\_i(\mathbf{x}\_i) - b \right] + \mu \left[ \sum\_{i=1}^{n} h\_i(\mathbf{x}\_i) - d \right]$$

$$\sum\_{i=1}^{n} g\_i(\mathbf{x}\_i) = b \tag{4}$$

$$\sum\_{i=1}^{n} h\_i(\mathbf{x}\_i) \ge d$$

where L represents the large system Lagrange function; Li represents the Lagrange function of the i-th subsystem; n represents the number of subsystems; xi represents the decision vector of the i-th subsystem (in this study, it represents water released from reservoirs); fi(xi) represents the evaluation function of the i-th subsystem; *n i*=1 *gi*(*xi*) = *b* is the equality constraint between subsystems; *n i*=1 *hi*(*xi*) ≥ *d* is the inequality constraint between subsystems; λ is the equality Lagrange multiplier constraint; and μ is the inequality Kuhn-Tacker multiplier constraint.

As shown in Equation (4), the equality constraint and inequality constraint are frequently used for constraining water balance and constraining cooperative work, respectively.

In this work, the subsystems' optimal solutions were solved using an adaptive genetic algorithm (AGA). The coordination between subsystems was implemented by means of coordinating variable iterations (Lagrange multiplier or Kuhn-Tacker multiplier) within the coordination level. The mixing method, which includes the interaction balance method and interaction prediction method, was employed for solving coordinating variables.

#### 2.2.2. Coordination between Multi-Scale Subsystems

It is clear that the subsystems' hydraulic connections and cooperative work between subsystems are often regarded as interconnection constraints between subsystems of the same scale. The multi-scale coupled ecological dispatching model accounts for the coordination between different scale subsystems. At first, the water balance constraint is the only interconnection constraint in the multiple hierarchical control system:

$$\sum\_{t=t\_0}^{t\_f} \mathbf{x}\_s(t) = \mathbf{x}\_i \tag{5}$$

where t0 and t*<sup>f</sup>* are the start and finish time of subordinate subsystems, respectively (smaller scale); xs represents the subordinate subsystem decision variable; and xi is the result of superior subsystem (larger scale).

As shown in Equation (5), the sum of the discharge flows in the small-scale model is equal to the discharge flow in the large-scale model for the corresponding period. Clearly, the water balance constraint is an equality constraint. Therefore, a new Lagrange multiplier, λ*s*, was introduced to regulate the water balance constraint. Usually, the Lagrange equation L is used as the whole basin scheduling model objective function. Thus, the mathematical expression of the multiple hierarchical control system is as follows:

$$L = \sum\_{i=1}^{n} \left\{ f\_i(\mathbf{x}\_i) + \sum\_{t=t\_0}^{t\_f} f\_\mathbf{s}(\mathbf{x}\_s) + \lambda\_s \left[ \sum\_{t=t\_0}^{t\_f} \mathbf{x}\_s(t) - \mathbf{x}\_i \right] \right\} + \lambda \left[ \sum\_{i=1}^{n} g\_i(\mathbf{x}\_i) - b \right] + \mu \left[ \sum\_{i=1}^{n} h\_i(\mathbf{x}\_i) - d \right]$$

$$\begin{array}{c} \sum\_{i=1}^{n} g\_i(\mathbf{x}\_i) = b \\ \sum\_{i=1}^{n} h\_i(\mathbf{x}\_i) \ge d \\ \sum\_{t\_f} \mathbf{x}\_f(t) = \mathbf{x}\_i \end{array} \tag{6}$$

where λ*<sup>s</sup>* is the Lagrange multiplier that regulates the water balance constraint between subsystems of different scales; *fs*(*xs*) represents the subordinate subsystem's (smaller scale) evaluation function belonging to the i-th subsystem. For our purposes, b represents the vector of the initial water amount in the reservoir at each time; d represents the vector of common water supply mission between reservoirs. The other symbols are the same as those defined in Equations (4) and (5).

As demonstrated in the above formulas, subsystems of different scales have different evaluation functions. However, solving a hierarchical control model with different objective functions for different layers is very difficult. As such, in this work, large-scale subsystem objective functions were converted to the subordinate subsystem's constraint conditions. In this way, objective functions of the large-scale subsystems were implemented by the subordinate subsystems. Unlike common multi-scale models that use results of the superior subsystem merely for subordinate subsystem boundary conditions, multi-scale subsystems—which are based on the theory of large-scale systems—consider coordination between different scale subsystems via an iteration of the coordination variables (Lagrange multiplier or Kuhn-Tacker multiplier).

The implication of coordination variables is deduced by the following interaction prediction method:

$$\text{Let}: \frac{\partial L}{\partial \lambda\_s} = \sum\_{t=t\_0}^{t\_f} \mathbf{x}\_s(t) - \mathbf{x}\_i = 0 \tag{7}$$

In order to satisfy Karush–Kuhn–Tucker conditions [35]:

$$
\frac{\partial \mathcal{L}}{\partial \mathbf{x}\_i} = \frac{\partial f\_i}{\partial \mathbf{x}\_i} - \lambda\_\mathbf{s} + \lambda \frac{\partial \mathcal{g}\_i}{\partial \mathbf{x}\_i} + \mu \frac{\partial h\_i}{\partial \mathbf{x}\_i} = 0 \tag{8}
$$

$$\frac{\partial L}{\partial \left[\sum\_{t=t\_0}^{t\_f} \mathbf{x}\_s(t)\right]} = \frac{\partial \left[\sum\_{t=t\_0}^{t\_f} f\_s(\mathbf{x}\_s)\right]}{\partial \left[\sum\_{t=t\_0}^{t\_f} \mathbf{x}\_s(t)\right]} + \lambda\_s = 0 \tag{9}$$

Thus:

$$\lambda\_s = \frac{\partial f\_i}{\partial \mathbf{x}\_i} + \lambda \frac{\partial \mathbf{g}\_i}{\partial \mathbf{x}\_i} + \mu \frac{\partial l\_i}{\partial \mathbf{x}\_i} = -\frac{\partial \left[\sum\_{t=t\_0}^{t\_f} f\_s(\mathbf{x}\_s)\right]}{\partial \mathbf{x}\_i} \tag{10}$$

As shown in the above formula, the large system global optimal solution is obtained exclusively at the moment when the derivative of the subordinate subsystem to the decision variable and the derivative of the superior subsystem to the decision variable are equal.

#### 2.2.3. Construction of the Subsystems

After a coordination level was assigned for the coordination variables, the subsystem was converted into a relatively simple problem. Hence, AGA was used to mathematically solve the subsystem. The subsystem's objective function is written as follows:

$$L\_i = f\_i(\mathbf{x}\_i) - \alpha \lambda\_i \mathbf{x}\_i + \lambda\_i g\_i(\mathbf{x}\_i) + \mu\_i h\_i(\mathbf{x}\_i) \tag{11}$$

where Li represents the Lagrange function of the i-th subsystem; *fi*(*xi*) represents the evaluation function of the i-th subsystem; α is a 0-1 variable, which is used as a subordinate subsystem structural discrimination coefficient, 1 and 0 denote whether the i-th subsystem has a subordinate subsystem or not, respectively; λ*<sup>i</sup>* and μ*<sup>i</sup>* refer to the corresponding coordination variable components of the i-th subsystem; and the other symbols are the same as those defined in Equations (4) and (5).

The mathematical expression of *fi*(*xi*) is as follows:

$$f\_i(\mathbf{x}\_i) = (1 - \alpha) \left( \sum\_{t=1}^{T} \frac{\mathbf{W} - \mathbf{x}\_i}{\mathbf{W}} + \sum\_{j=1}^{M} D\_j \right) + \alpha \sum\_{j=1}^{N} D\_j \tag{12}$$

where W refers to the subsystem water demand; xi refers to the subsystem water supply; *Dj* refers to the degree of change of the j-th index; M and N represent the number of superior subsystem and subordinate subsystem indicators, respectively; and the other symbols are the same as those defined in Equations (4) and (5).

#### *2.3. A Solution of the Multi-Scale Coupled Ecological Dispatching Model*

The gradient method was used to mathematically resolve the coordination level. The coordination variables were all initialized to zero; and the subsystem's calculation results were passed to the coordination level. The iterative formulas are as follows:

$$\sum\_{i=1}^{n} g\_i(x\_i) = b$$

$$\sum\_{t=t\_0}^{t\_f} x\_i(t) = x\_i \tag{14}$$

$$
\lambda\_s^{j+1} = -\frac{\partial \left[ \sum\_{t=t\_0}^{t\_f} f\_s(\mathbf{x}\_s) \right]}{\partial \mathbf{x}\_i} \tag{15}
$$

$$\mu\_i \mu^{j+1} = \mu\_i + [h\_i(\mathbf{x}\_i) - d] \ast \mathbf{R} \tag{16}$$

$$
\lambda\_i \lambda\_i^{j+1} = \frac{\partial f\_{i+1}}{\partial \mathbf{x}\_i} - \lambda\_{\mathbf{s}\_{i+1}} + \lambda\_{i+1} + \mu\_{i+1} \tag{17}
$$

where the symbol superscripts represent the number of iterations; unmarked symbols represent symbols of the j-th iteration; symbol subscripts represent the number of subsystems; R refers to the account step of the coordination variables; <sup>∂</sup> *fi*+<sup>1</sup> <sup>∂</sup>*xi* represents the derivative of the objective function of the i + 1-th subsystem to initial water of this subsystem (equal to surplus water of the last subsystem); ∂[ *t* - *f fs*(*xs*)]

*t*=*t* 0 <sup>∂</sup>*xi* represents the derivative of the objective function of the inferior subsystem to the superior subsystem output; and the other symbols are the same as those defined in Equations (4)and (5).

After the calculation, the coordination variables were passed to the subsystems until the convergence condition was met. The coordination variables' convergence conditions are expressed as:

$$\left|\lambda\_s^{j+1} - \lambda\_s^{j}\right| \le \varepsilon\_{\lambda\_s} \tag{18}$$

$$\left|\mu\_{l}^{j+1} - \mu\_{l}^{j}\right| \leq \varepsilon\_{\mu\_{i}} \tag{19}$$

$$\left|\lambda\_i^{j+1} - \lambda\_i^{j}\right| \le \varepsilon\_{\lambda\_i} \tag{20}$$

where ελ*s*, εμ*<sup>i</sup>* , ελ*<sup>i</sup>* refer to the coordination variables' convergence thresholds.

Unlike the traditional targets, the derivative of the ecological objective function to initial water is difficult to calculate. Thus, in this work, a differential response model between the ecological objective function and initial water of a small-scale subsystem was established. The results were calculated in advance for coordination level.

#### *2.4. Description of the Study Area*

The lower reaches of the Yellow River (Figure 1) refer to the Yellow River stream segment below Peach Blossom Valley. At this location, the river length is 786 km and the drainage area covers 23,000 km2. Furthermore, the stream gradient is small and the river's ecological problems are serious. The main water conservancy projects within this area are the Xiaolangdi reservoir and the Sanmenxia reservoir; which effectively store 5.1 billion cubic meters (bcm) and 0.462 bcm of water, respectively. Annual consumptive water use in the study area is 7.5 bcm. During flood season, there are water resources to ensure the implementation of ecological operation.

**Figure 1.** Locations of the downstream region of the Yellow River.

The ecological problems in the lower reaches of the Yellow River differ depending on the time of year. (1) During the major flood period—from July to September—the main ecological problem is the disappearance of the flood and high flow pulses, resulting in a flat, non-fluctuating flow; (2) during the drought period—from November to March—the river suffers a water shortage and fails to meet the minimum instreaming ecological flow requirements; (3) during medium season, the primary problem is habitat change in the riverine biological communities.

Because the characteristic flow regimes are relatively complete in the major flood period, for this study, the major flood period was selected as the ecological critical period that was assessed using the small-scale scheduling model. Water demand data used in this work came from the Comprehensive Planning of the Yellow River Basin [36], which is a book compiled by the Yellow River Conservancy Commission (YRCC) of the ministry of water resources of China; and the daily runoff data was taken from the Hydrological Yearbook of the People s Republic of China [37].

The multi-scale coupled ecological dispatching model shown in Figure 2 was established based on the characteristics of the study area described above.

In Figure 2, x refers to the discharge from the reservoirs; λ<sup>1</sup> and λ<sup>2</sup> are the Lagrange multipliers that regulate the water balance constraint between the Xiaolangdi and Sanmenxia reservoir superior subsystems, respectively, at different times; λ<sup>3</sup> and λ<sup>4</sup> are the Lagrange multipliers that regulate the water balance constraint between the superior subsystems and subordinate subsystems, respectively, of the Xiaolangdi and Sanmenxia reservoirs; and μ is the Lagrange multiplier that regulates the water balance constraint between the superior subsystems during the same period.

The objective function of the model is as follows:

(1) Minimum water shortage rate:

$$\text{Minf}\_{S-D} = \sum\_{t=1}^{T} \sum\_{m=1}^{M} \frac{\mathbf{D}\_{(m,t)} - \mathbf{S}\_{D(m,t)}}{\mathbf{D}\_{(m,t)}} \tag{21}$$

where f*S*−*<sup>D</sup>* refers to the objective function of the water deficient ratio; m refers to the serial number of the reservoir; M refers to the sum of reservoirs; t refers to the scheduling period number; T refers to scheduling period length; D(*m*,*t*) refers to water demand of reservoir m at the t-th period; and S*D*(*m*,*t*) refers to the water supply of D(*m*,*t*);

#### (2) The minimum degree of hydrologic alteration is expressed as:

$$\text{Min}\sum\_{i=1}^{I} D\_i = \sum\_{L=1,2,3} \left| \frac{N\_i(La) - N\_c(La)}{N\_c(La)} \right| \tag{22}$$

where I represents the total number of indicators; La = 1, 2, 3 represents the high, middle and low target ranges, respectively; *Di* refers to the degree of change of the i-th index; *Ni*(*La*) refers to the number of years in which the i-th index is observed in the corresponding target ranges L; *Ne*(*La*) refers to the number of years in which the i-th index is expected to be observed in the corresponding target ranges L.

The objective function (Lagrange function) of the overall system model is obtained by substituting Formula (21) and Formula (22) into Formula (6).

Constraint conditions of the model include water level constraint, flow constraint, water balance constraint, and flow balance constraint.

(1) Water level constraint:

$$Z\_{\min}(m, t) \le Z(m, t) \le Z\_{\max}(m, t) \tag{23}$$

(2) Flow constraint:

$$Q\_{\rm Omin}(m, t) \le Q\_{\rm O}(m, t) \le Q\_{\rm Omax}(m, t) \tag{24}$$

(3) Water balance constraint:

$$V(m, t+1) = V(m, t) + Q\_l(m, t) - Q\_\bullet(m, t) \tag{25}$$

(4) Flow balance constraint:

$$Q\_l(m+1,t) = Q\_\mathcal{O}(m,t) + q(m,t) \tag{26}$$

(5) Water supply reliability constraint:

$$\begin{array}{l}reliability\_{algorithm} \ge 0.75\\reliability\_{urban} \ge 0.95 \end{array} \tag{27}$$

where *Z*(*m*, *t*) refers to the water level of reservoir m at period t; *Z*max(*m*, *t*), *Z*min(*m*, *t*) are the top and bottom limitations of the water level of reservoir m, respectively at period t; *QO*(*m*, *t*) refers to the discharged flow of reservoir m at period t; *QO*max(*m*, *t*), *QO*min(*m*, *t*) are the top and bottom limitations of the discharged flow of reservoir m at period t, respectively; *V*(*m*,*t*) refers to the storage of reservoir m at period t; *QI*(*m*,*t*) refers to the inflow of reservoir m at period t; *q*(*m*,*t*) refers to the local inflow of reservoir m at period t; and *reliabilityagricultural* and *reliabilityurban* refer to the reliability of agricultural and urban water supply, respectively (urban water supply includes the industrial water supply and domestic water supply). The reliability data comes from the Comprehensive Planning of the Yellow River Basin.

Among the constraint conditions, water balance constraint between scheduling periods and flow balance constraint were implemented by coordination variables.

In this study, the daily flow from 1958 to 1977 was used as the reference series and the daily flow from 2007 to 2014 was selected as the scheduling period. The reservoir inflow, domestic water demand, agricultural water demand, and ecological water demand, etc., are required inputs of the multi-scale coupled ecological dispatching model.

**Figure 2.** Information transfer structure diagram of the model.

#### **3. Results**

For comparison purposes, the hydrologic alterations caused by the operation of the current reservoir in the middle and lower Yellow River were assessed using the RVA method. The daily flow from 1958 to 1977 was used as the pre-impact series, and 2007–2014 were selected as the post-impact series. The degree to which the RVA target range is not attained is a measure of hydrologic alteration and is calculated as follows:

$$\langle | \text{Observed frequency} - \text{Expected frequency} \rangle \langle \text{Expected frequency} \tag{28}$$

The RVA uses 33 hydrological parameters to evaluate the hydrologic alterations [4], which are categorized into five groups: magnitude of monthly water conditions, magnitude and duration of annual extreme conditions, timing of annual extreme water conditions, frequency and duration of high and low pulses, and rate and frequency of water condition changes. As shown in Figure 3, the integrated hydrologic alteration reaches 86%. In terms of indicators, most change was observed in the magnitude of monthly water conditions and the magnitude and duration of annual extreme conditions. Among those, the change in magnitude of monthly water conditions during the flood season is the largest. In terms of RVA categories, the high RVA categories reduced the most.

**Figure 3.** Hydrologic alteration of current scheduling.

The multi-scale coupled ecological dispatching model was established to generate small-scale ecological flow regimes and avoid dimension disaster. In this study, subsystems were solved by AGA and coordinated via the coordination variables. As such, the degree of hydrologic alteration can be used as an evaluation function of the reservoir scheduling model, which is unattainable for the conventional models. Thus, the multi-scale coupled ecological dispatching model provides more information to reservoir managers. Compliments of decomposition-coordination algorithms and AGA, large-scale model and small-scale model results can efficiently converge to the optimal solutions within 200 iterations as shown in Figure 4.

**Figure 4.** (**a**) Large-scale model optimization processes; (**b**) Small-scale model optimization processes.

Figure 5 compares the simulative monthly flow and actual monthly flow from 2007 to 2014. The results show that the model is capable of solving the challenges associated with uneven spatial and temporal distribution of water resources. Three observations are evident in Figure 6: (1) The simulative monthly flow is more fluctuant than the actual monthly flow; (2) the simulative monthly flow in the flood season is larger than the actual monthly flow; and (3) the occurrence time of the simulative high flow is more regular than actual high flow. These results are explicable, as more water was needed during the ecologic critical period to generate the characteristic flow regimes. Moreover, the average flow in flood season can be increased to 1400 m3/s under the premise of meeting the water supply reliability.

**Figure 5.** Monthly discharge from Xiaolangdi reservoir.

The 'percent-of-flow' (POF) approach has been widely used in the USA, the European Union, and elsewhere over the past decade [3]. To contrast the above results, the flow was also generated using the POF approach, which sets protection standards by using allowable departures from natural conditions, expressed as percentage alteration. The hydrologic alteration of flow was calculated under the POF approach (Figure 6) and contrasted with the hydrologic alteration of simulative flow (Figure 6). Comparison of Figures 3 and 6 shows that the POF approach is more protective with respect to the magnitude of monthly water conditions—as the Hydrologic alteration of the magnitude of monthly water conditions were reduced from 81% to 65%. Yet, in other groups, particularly in the magnitude and duration of annual extreme conditions, the POF approach does not provide a high degree of protection for natural flow variability. Contrarily, the hydrologic alteration of the simulative flow was reduced in all five groups. The integrated hydrologic alteration was reduced from 86% to 53%. On the one hand, this is because the objective function of the optimization model is to reduce the degree of hydrological alterations. On the other hand, the incoming water in the whole dispatching period of the optimization model is known, and this is not possible in the actual scheduling process. The hydrologic alterations of five groups are 0.60, 0.41, 0.29, 0.67, and 0.67, respectively (see the beginning of this section for corresponding groups). It is worth noting that the simulative flow did not provide a very high degree of protection for the frequency and duration of high pulses, particularly in the high RVA categories. This is explained by the difficulties incurred with providing protection for natural flow variability of different indicators in different RVA categories. For instance, more water was needed to keep the frequency of flood pulses in the high RVA category; and the increased water concentration reduced the unit water value, resulting in a lack of water elsewhere.

From the perspective of hydraulic parameters, the number of days within the appropriate velocity range during the reference series (1958–1977) is 2117, which is 29% of the total number of days. The number of days within the appropriate velocity range during the scheduling period (2007–2014) is 1452, which is 50% of the total number of days. Due to the influence of flow flattening, the number of days suitable for fish has increased. We expected that the duration of the simulated flow within the appropriate velocity is not less than the duration of natural flow. The number of days within the appropriate velocity range in the simulated flow accounted for 46% of the total. And the number of days within the appropriate velocity range under the POF approach accounted for 43% of the total. Obviously, this is an easy goal to achieve.

To solve the problem of uncertainty, we use the autoregressive model to generate nine sets of inflow runoff series. Ten sets of inflow runoff series (Including the actual inflow series from 2007–2014) were used as input for simulation calculation. The final result is the average of 10 sets of calculated results and the variation range of 10 sets of results. The values above are the mean, and the values in

brackets are the variation ranges. The quantitative evaluation and biological significance of various indicators under the POF approach and optimization model respectively are shown in Table 3. As it is difficult to simulate decision-makers' decisions, the uncertainty of current scheduling is not shown. The hydrological alterations of the current scheduling are only used as a reference in Table 3. We found that the calculation results of the optimization model performed well in the second, third and fifth categories of indicators. That means the model performed better in terms of vegetation expansion, fluvial topography, natural habitat construction, fish migration, and spawning and life cycle reproduction compared to current scheduling and the POF approach. In terms of uncertainty, the optimization model performs better than the POF approach too. Especially in the second, third and fifth categories of indicators, the performance of the optimization model are far better than the POF approach, and the variation ranges barely overlap. That means the model performed better in addressing uncertain problems compared to current scheduling and the POF approach in terms of the above ecological problems.

**Figure 6.** (**a**) Hydrologic alteration of POF approach; (**b**) Hydrologic alteration of simulative flow.



Table 3 shows the quantitative evaluation and biological significance of various indicators under the POF approach and optimization model. The values of the first five categories indicators are degrees of hydrologic alterations. The value of the last indicator represents the proportion of the number of days in the suitable ecological velocity range to the total number of days.

#### **4. Discussion**

Compared to conventional models, the multi-scale coupled ecological dispatching model can obtain the degree of hydrologic alteration of small-scale ecological flow regimes and daily stream flow. At the same time, the results can be manipulated by adjusting the weights afforded to the degree of hydrologic alteration and other targets. In contrast, degrees of hydrologic alteration based on whether or not to consider the characteristic flow requirements and whether or not to introduce the small-scale model are shown in Table 4. When considering the characteristic flow requirements, reservoirs are required to release more water in the ecological critical period. Under these circumstances, the appropriate characteristic flows can be generated by the small-scale model. When the characteristic flow requirements are not taken into consideration, the agricultural water supply dependability increases from 75% to 90%.


**Table 4.** Contrast of the degree of hydrologic alteration.

Table 4 shows the influence of considering characteristic flow requirements on guaranteed water supply reliability and hydrologic alteration.

Compared with the existing research, the existing model either directly optimizes the scheduling rules [17] or evaluates the existing daily flow process [24]. However, the rule optimization model cannot describe the small scale runoff process in detail, and the runoff evaluation model cannot directly guide reservoir operation. Therefore, a modelling approach which allows scheduling of reservoir operation in such a way that hydrological alteration is reduced was presented in this work. The multi-scale coupled ecological dispatching model provides a framework for the formulation of multi-scale scheduling rules. The more flexible decoupling method can be used by the model to optimize the scheduling rules directly. In further research, the large-scale model can optimize the operation rules, while the small-scale model can optimize the discharge process of reservoirs. This method is suitable for areas where the hydrological regime is seriously affected by dam operation. This model is more suitable for large rivers because of its attention to the degree of hydrological alterations. For other river basins, decision makers can choose the key ecological issues of the river basin, and the weight of the objective function can be formulated according to the key ecological problems and the ecological significance of various indicators. Scheduling rules suitable for this basin can be found by this way.

However, the multi-scale coupled ecological dispatching model established in this study still has some issues. The competitive relationship between the ecological and economical target was coordinated by the weighting method. A more reasonable and flexible multi-objective coordination mechanism should be introduced to the model. Furthermore, considering the impact of hydrological uncertainty, calculation results from the model cannot directly guide the reservoir operation. A multi-scale coupled scheduling rule for the reservoir should be put forward and planned for future studies, which is what we are going to do next.

In general, the multi-scale coupled ecological dispatching model is capable of considering the degree of hydrologic alteration of small-scale ecological flow regimes. The ecological target was regarded as both a water requirement and a difference between the artificial flow and natural flow. Operation rules extracted from results of the coupled ecological dispatching model will include more information and be more reasonable.

#### **5. Conclusions**

A multi-scale coupled ecological dispatching model, that considers the degree of hydrologic alteration, was established to satisfy the multi-scale features of ecological water demand. In order to acquire all characteristics of the environmental flows and to avoid the "dimension disaster", the decomposition coordination algorithm was applied to this model. Furthermore, this model is capable of using account steps from months to hours, as needed. Moreover, the general multi-scale model uses results of the superior subsystems only as boundary conditions of the subordinate subsystems; while in our model, multi-scale subsystems are based on the theory of large-scale systems and we consider the coordination of multi-scale objectives.

Simulation results show that the model is able to solve the problem of uneven spatial and temporal distribution of water resources and obtain degrees of hydrologic alteration of characteristic flow regimes. The difficulty of generating characteristic flow regimes is indicated by the results, which will be helpful to the formulation of operation rules. The calculation results indicate that increasing the average flow in flood season to 1400 m3/s can reduce the degree of hydrologic alteration from 86% to 53% in the downstream Yellow River. Model optimization results were compared with results of the POF approach. The comparison showed that further reduction in hydrologic alteration is possible, but further research is needed to determine which method is more suitable for accomplishing this task.

While the proposed model is overall a success, there are some limitations worth noting. On the model side, a more reasonable and flexible multi-objective coordination mechanism needs to be established. From the mathematics side, a high-efficiency and applicative arithmetic method needs to be developed or improved to extract scheduling rules from artificial daily flow series.

**Author Contributions:** Conceptualization, T.Z. and Z.D.; methodology, T.Z.; software, T.Z.; investigation, X.G. and Z.H.; resources, X.G. and R.S.; data curation, R.S.; writing original draft preparation, T.Z.; writing—review and editing, T.Z.; W.W. and Z.D.

**Funding:** This research has been financially supported by the National Key Research and Development Program of China (2018YFC1508200), the Postgraduate Research & Practice Innovation Program of Jiangsu Province (KYCX18—0574) and the Fundamental Research Funds for the Central Universities (2018B603X14).

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **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* **SPH Simulation of Interior and Exterior Flow Field Characteristics of Porous Media**

#### **Shijie Wu 1, Matteo Rubinato <sup>2</sup> and Qinqin Gui 1,\***


Received: 16 January 2020; Accepted: 17 March 2020; Published: 24 March 2020

**Abstract:** At the present time, one of the most relevant challenges in marine and ocean engineering and practice is the development of a mathematical modeling that can accurately replicate the interaction of water waves with porous coastal structures. Over the last 60 years, multiple techniques and solutions have been identified, from linearized solutions based on wave theories and constant friction coefficients to very sophisticated Eulerian or Lagrangian solvers of the Navier-Stokes (NS) equations. In order to explore the flow field interior and exterior of the porous media under different working conditions, the Smooth Particle Hydrodynamics (SPH) numerical simulation method was used to simulate the flow distribution inside and outside a porous media applied to interact with the wave propagation. The flow behavior is described avoiding Euler's description of the interface problem between the Euler mesh and the material selected. Considering the velocity boundary conditions and the cyclical circulation boundary conditions at the junction of the porous media and the water flow, the SPH numerical simulation is used to analyze the flow field characteristics, as well as the longitudinal and vertical velocity distribution of the back vortex flow field and the law of eddy current motion. This study provides innovative insights on the mathematical modelling of the interaction between porous structures and flow propagation. Furthermore, there is a good agreement (within 10%) between the numerical results and the experimental ones collected for scenarios with porosity of 0.349 and 0.475, demonstrating that SPH can simulate the flow patterns of the porous media, the flow through the inner and outer areas of the porous media, and the flow field of the back vortex region. Results obtained and the new mathematical approach used can help to effectively simulate with high-precision the changes along the water depth, for a better design of marine and ocean engineering solutions adopted to protect coastal areas.

**Keywords:** Smooth Particle Hydrodynamics (SPH); porous media; mathematical model; coastal structure; ocean and engineering

#### **1. Introduction**

In ocean and water conservancy projects, pore structures are widely used, such as reservoir dams, artificial reefs, and breakwaters. When water flows through pore structures, factors such as porosity rate and volume of porous media can be the major cause for different flow field fluctuations [1]. Therefore, to investigate the influence of porosity rate and volume of porous media on the flow field, which is crucial to characterize the features of flow field characteristics around the porous media [2], this paper analyzes the internal and external flow field characteristics of typical structures by simulating these physical phenomena.

The commonly used simulation method to study these effects is CFD (Computational Fluid Dynamics). CFD simulation can effectively simulate the flow field characteristics, however multiple complications and lengthy simulations are associated with the establishment of a proper mesh [3,4]. Boundary element simulations can effectively simulate the influence of water flow on pore structure, but these methods only divide the elements on the boundary domain to meet the approximations applied to the functions implemented for the boundary conditions of the control equation. Therefore, compared with the finite element method, the boundary element method has the advantage of fewer elements to consider and a simpler data preparation. However, when the boundary element method is used to solve the nonlinear problem, which corresponds to the integral corresponding to the nonlinear term, more complications are frequently encountered because this kind of integration has strong singularity near single points, making an accurate solution very challenging to be obtained. Therefore, to reduce and avoid these complications, the key properties of the Smooth Particle Hydrodynamics (SPH) method, which is based on the pure Lagrangian numerical simulation method, were adopted for this study. This method can track the trajectory of each particle, and has good advantages in terms of visualization, high precision, and large processing. Many studies have been completed to date to explore the hydrodynamic characteristics of porous media. For example, Hossei Basser et al. [5] used SPH method to characterize the gravity flow of a single fluid in porous media and multi-fluid flows in impermeable layers and porous media, as well as the influence of unidirectional flow on pore structure and particle motion. The results that were obtained confirmed how the SPH numerical simulation method had the characteristics of high accuracy and visualization, which could effectively and accurately simulate the movement of each particle in the flow field inside and outside the pore structure, describing very accurately the characteristics of the flow field inside and outside the pore structure.

Pore flow field simulation can be simulated focusing on two aspects: internal flow field and external flow field, which refer to the specific flow field formed in large-scale ocean engineering, such as artificial reefs, dam break and breakwater, under the condition of laminar and turbulent flows. Boundary conditions are very important and are typically periodic or open. The first category means that the flow particles run from the inlet boundary to the outlet boundary, re-entering the inlet boundary and form a continuous loop. The second category refers to the flow particles moving in from the inlet boundary to the outlet boundary without re-entering from the inlet boundary, because this approximation considers that new flow particles will move in following an uninterrupted process. Under the open boundary conditions, researchers typically simulate the internal and external flow fields, effectively realizing the continuous process of particles motion [6]. However, by using these conditions although the accuracy of the inlet and outlet flows is high, the interaction with the porous media is relatively complex and affected by lower accuracy. Therefore, periodic boundary conditions accompanied with a damping zone were employed in this paper to complete the work mentioned above.

When simulating pore structure models, porosity must be considered, because factors such as pore permeability can affect the flow field around the pore, like for example in processes such as dam break, where the dam body is composed of pore structure. Therefore, under different porosity, the dam break strength is different, generating different impacts in case of failure of the dam taken into consideration. Additionally, when the flow comes in contact with the pore structure, the seepage characteristics of the pore structure can further affect the surrounding flow field [7], so the porosity that affects the pore structure is one of the key factors to be analyzed.

To date, multiple studies have been conducted on the characterization of pore structures. Qinqin Gui et al. [8] studied the coupling law of wave and pore structure by using the incompressible SPH method. Through numerical simulation and experiment, the typical flow around and inside the underwater porous structure was compared, however the analysis of the internal and external flow field of the porous structure was not considered as well as the influence of water flow on the porous structure. Abbas Khayyer et al. [9] proposed a numerical solution for variable porosity, because the purpose of this study was to investigate the effect of waves on structures with different porosity. Vanneste and Troch [10] conducted a detailed two-dimension numerical study on the wave propagation and verified the calculation results in terms of average water level and spatial pore pressure distribution from the incident wave field to the normal wave field. Shao [11] proposed an Incompressible Smooth Particle Hydrodynamics (ISPH) method to simulate the interaction between waves and porous media,

and achieved interesting results on the analysis of flows and pore structures: through comparative analysis, the error obtained by Shao [11] was between 10% and 15%, which better reflects the accuracy of numerical simulations when the SPH is applied [8]. To summarize what was found in literature, to date, traditional SPH method is generally used to estimate the influence of waves on the porous media, but seldom for the coupling the effects of water flow and pore separately. Therefore, there is a need for research to be conducted to demonstrate that the SPH method can be adopted to simulate the action process of water flow on porous materials and the consequent effects. To fulfill this gap, SPH is used in this work to simulate coupling effects of water flow and porous media with the periodic circulation flow as a boundary condition. The influence of the flow field on the porous media was analyzed, together with the flow field inside and outside the porous media, the longitudinal and vertical velocity distribution, and the eddy current law of the eddy current field.

#### **2. Smooth Particle Hydrodynamics (SPH) Methodology**

The Smooth Particle Hydrodynamics method [12–22] is commonly used to describe a continuous fluid (or solid) as a group of interacting particles, despite each of these particles carries independent physical features such as different mass and dissimilar velocity. This method tends to characterize the continuous fluid by solving the behavior of the entire group of particles to simplify the problem and to determine the mechanical behavior of the entire system. In principle, as long as the number of particles is sufficient, the mechanical process can be accurately described. In practice, these particles are distributed irregularly, but in theory, they all move according to the motion law of the conservation governing equation. Any field of particles can be represented by applying a sufficiently accurate integral obtained by approximating the interpolation between adjacent particles. Considering this aspect, spatial derivatives such as gradient and scatter operators can be similarly estimated using the properties of all the particles in the Navier-Stokes equations, and the kernel approximation is commonly used to approximate the field function.

The procedure applied for this research uses the Smooth Particle Fluid Dynamics (SPH) method, which can be approximated by integral interpolation:

$$A(r) = \int A(r')\mathcal{W}(r - r', h)dr'\tag{1}$$

where *A*(*r*) represent any function. For numerical simulations, the integral interpolation is calculated as follows:

$$A\_{a}(r) = \sum\_{b} {}\_{b} m\_{b} \frac{A\_{b}(r)}{\rho\_{b}} {}\_{ab} \tag{2}$$

where *a* and *b* represent particles; *mb* and ρ*b* represent particle mass and density, *Wab* is a kernel function, r is the distance between the reference particle and the adjacent (m), h is the smoothing length (m). To ensure a linear and angular momentum conservation, the asymmetric pressure gradient terms can be obtained as follows:

$$\left(\frac{1}{\rho}\nabla P\right)\_a = \sum\_b \, m\_b (\frac{P\_a}{P\_{a^2}} + \frac{P\_b}{P\_{b^2}}) \nabla\_a \mathcal{W}\_{ab} \tag{3}$$

where ∇*aWab* represents the kernel gradient taken relative to the position of the particle *a* and P represents the pressure; similarly, the divergence of the vector at a given particle *a* can be estimated by *u*:

$$
\nabla \mu\_a = \rho\_a \sum\_b \rho\_b m\_b (\frac{\mathcal{U}\_a}{P\_{a^2}} + \frac{\mathcal{U}\_b}{P\_{b^2}}) \nabla\_a \mathcal{W}\_{ab} \tag{4}
$$

The viscous force can be calculated by combining the gradient and the divergence term, as follows:

$$(\nu\_0 \nabla^2 u)\_a = \sum\_b \, m\_b \frac{2(\nu\_a + \nu\_b)}{(\rho\_a + \rho\_b)} \cdot \frac{\mu\_{ab} r\_{ab} \cdot \nabla\_a W\_{ab}}{|r\_{ab}|^2 + \eta^2} \tag{5}$$

where υ represents the viscosity coefficient. Since pressure is highly sensitive to turbulence, it has been necessary to adopt the Laplace term in the pressure Poisson equation, written as follows [23]:

$$\nabla \cdot \left(\frac{\nabla P}{\rho}\right)\_a = \sum\_b \, \_b m\_b \frac{8 \cdot P\_{ab} r\_{ab} \cdot \nabla\_a W\_{ab}}{\left(\rho\_a + \rho\_b\right)^2 \cdot \left(\left|r\_{ab}\right|^2 + \eta^2\right)}\tag{6}$$

$$P\_{ab} = P\_a - P\_b \ \ r\_{ab} = r\_a - r\_b \tag{7}$$

The kernel function uses the cubic spline kernel function developed by Monaghan [24], which is related to the dimension of the field investigated and the distance between the relevant points. The kernel function is displayed below:

$$\mathcal{W}(r,h) = \begin{cases} \left(10 \times \left(1 - 1.5q^2 + 0.75q^3\right)\right) / 7\pi h^2 & 0 \le q < 1\\ \left(10(2 - q)^3\right) / 28\pi h^2 & 1 \le q \le 2\\ 0 & q > 2 \end{cases} \tag{8}$$

where *r* is the distance between the reference particle and the adjacent one; *q* = *r*/*h* is the relative distance.

#### **3. SPH Model**

#### *3.1. Equations for Flow Field Porous Media*

Based on the previous studies conducted [19,25], the flow outside of the porous media is typically considered laminar and can be solved by the two-dimensional Navier-Stokes equations. The external flow field of pore structure is as follows:

$$\frac{d\rho}{\rho dt} + \nabla \cdot \boldsymbol{u}\_{\rm w} = 0 \tag{9}$$

$$\frac{du\_w}{dt} = -\frac{1}{\rho}\nabla \mathbf{P} + \mathbf{g} + \ \nu \nabla^2 u\_w \tag{10}$$

where ρ is the flow particle density; *t* is the time; *uw* is the external flow rate of the pore.

Equation (10) is the governing equation of the external flow field of the porous media. This equation incorporates the influence of the current on the porous media during the movement, and the parameters to calculate the motion of the eddy current.

The internal flow field inside pore structure can be obtained as follows [19,25]:

$$\frac{du\_p}{dt} = -\frac{1}{\rho}\nabla P + \mathbf{g} + \boldsymbol{\nu} \cdot \nabla^2 u\_p - \frac{C\_f n\_w^2}{\sqrt{\mathbf{K}\_P}} u\_p |u\_p| \tag{11}$$

$$K\_p = 1.643 \times 10^{-7} \left[ \frac{d\_{50}}{d\_0} \right]^{1.57} \left[ \frac{n\_w^{-3}}{\left(1 - n\_w\right)^2} \right] \text{ C}\_f = 100 [d\_{50}(m) \cdot \left(\frac{n\_w}{K\_p}\right)^{1/2}] \tag{12}$$

where *up*<sup>1</sup> is the conveying speed (*m*/*s*); *nw* is the porosity (/); *Kp* is the permeability (/); *Cf* is the nonlinear resistance coefficient (/).

#### *3.2. Numerical Model Solving Process*

In the first prediction step, the velocity and the position of the particles are calculated according to the momentum equations, and these equations can be expressed as follows [11]:

$$
\Delta u\_\* = \left(\mathbf{g} + \nu \nabla^2 u\_p - \frac{\nu \mathbf{n}\_w}{K\_P} \mathbf{u} - \frac{\mathbb{C}\_f n\_w}{\sqrt{K\_P}} \mathbf{u}|\mathbf{u}|\right) \Delta t \tag{13}
$$

$$
\Delta u\_\* = u\_t + \Delta u\_\* \tag{14}
$$

$$r\_\* = r\_l + \mu\_\* \Delta t \tag{15}$$

where Δ*u* is the particle velocity that changes during the prediction step; Δ*t* is the time increment; *ut* and *rt* represents the particle velocity and position at time t, respectively. The pressure term is based on the classical pressure Poisson equation that can be expressed as follows [23]:

$$\nabla \cdot \left(\frac{1}{\rho} \nabla P\_{t+1}\right) = \frac{\rho\_0 - \rho\_\*}{\rho\_0 \Delta t^2} \tag{16}$$

where ρ<sup>0</sup> represents the initial constant particle density; ρ\* represents the central particle density after the prediction step, and P*t*+<sup>1</sup> is the pressure of the particles at the *t*+1 time step. In the calibration of the second step, the pressure gradient term is combined with the momentum equation to ensure incompressibility. Pressure can be used to correct particle velocity as follows:

$$
\Delta u\_{\ast \ast} = \frac{-1}{\rho\_{\ast}} \nabla p\_{t+1} \Delta t \tag{17}
$$

$$
\mu\_{t+1} = \mu\_\* + \Delta u\_{\*\*} \tag{18}
$$

$$r\_{t+1} = r\_t + \frac{u\_{t+1} + u\_t}{2} \Delta t \tag{19}$$

where *ut*+<sup>1</sup> is represents the particle velocity and *rt*+<sup>1</sup> represents position at the moment of *t* + 1.

#### *3.3. Boundary Conditions*

#### 3.3.1. Free Surface Boundary

Shao and Lo [23] proposed a free surface treatment method where if the density of a particle is 10% lower than the reference density, it can be considered as a surface particle, and then zero pressure can be applied as a known boundary condition [8].

These assumptions are adopted for this study.

#### 3.3.2. Fluid-Structure Coupling Boundary

The interface region is placed between the exchange flow and the porous flow region at the boundary of the pure fluid interface, and the interface line is located at the center of the region. Therefore, the width of the two regions (the fluid portion and the porous portion) at the interface boundary line is twice the distance between particles. By using the proposed SPH interface model, the average value of the interface region is obtained by solving the pressure Poisson equation in each calculation time step using the SPH kernel function. Only the adjacent particles located in the interface area are included in the summation.

At the same time, the interface boundary conditions can automatically satisfy the normal and tangential continuity.

#### 3.3.3. Impermeable/Fixed Solid Wall Boundary

The impermeable/fixed solid wall is maintained along the boundary line. The pressure gradient between particles is adopted and when fixed virtual particles are used, the anti-skid boundary is selected.

3.3.4. Periodic Inflow and Outflow Boundaries Accompanied with a Damping Zone

According to previous studies conducted on pore structure logistics field, this paper adopts the boundary of inflow and outflow as displayed in Figure 1 below.

**Figure 1.** Flow boundary diagram.

Particles leaving one side of the computational domain move directly through the other side, and particles near the open side boundary interact with particles near the complementary open side boundary. The fluid particles near the left end of the channel can also interact with the fluid particles near the right end in the core affected area. In this sense, the upstream and downstream boundaries can actually be considered as overlapping. In the SPH numerical algorithm, we used the particle linked list to search for neighborhoods through a square grid with a side length of 2\*H. Therefore, the fluid particles in the first column of the block located near the upstream boundary can be in contact with the particles in the last column of the block located near the downstream boundary, and vice versa, so the exit boundary flow velocity is also consistent with the inlet boundary flow velocity trend. As a preliminary test to verify this hypothesis, we set the calculation area as 3 m, the inlet flow rate as 0.00936 m3/s, and the particle size D as 0.1 m. In order to further test the performance of our method, three time stages t = 18.0 s, t = 20.0 s, and t = 22.0 s were considered to analyze its inlet and outlet flow rate, respectively. x = 0.0 represents the inlet flow rate distribution curve; x = 60.0 represents the outlet flow rate distribution curve. Analyzing the tests shown in Figure 2, it is possible to notice that the flow distribution of the particles is very similar between inlet (x = 0.0) and outlet (x = 60.0). Furthermore, the velocity ranges from 0 *m*/*s* to 1.5 *m*/*s* for each time step displayed proving the periodic boundary conditions of the inlet and outlet.

**Figure 2.** Comparison of flow velocities at inlet and outlet boundary at different times at T = 18.0 s, T = 20.0 s and T = 22.0 s.

It can also be observed from Figure 2 that the inlet and the exit boundary particles are basically consistent in the horizontal direction despite different time steps considered, confirming that the particles at the exit layer can smoothly enter the inlet across the boundary.

However, if the flow is significantly disturbed before exiting the outlet boundary, the inflow condition will inherit these disturbances and affect the downstream flow. In that case, at least a damping zone is needed downstream the inlet boundary and a numerical calibration is needed to

clarify the effectiveness of the damping zone [26,27]. For this study, a damping zone is set at the junction of the exit boundary and the entrance boundary.

Figure 3 displays the distribution of the velocity of the flow at the exit boundary along the depth of the water after passing through the damping zone. It can be seen that the speed variations remain relatively stable along the depth of the water at different times (t = 60 s; t = 70 s; t = 80 s; t = 90 s), hence the flow velocity at the outlet remains basically identical (the average flow velocity is stable at about 0.13 *m*/*s* at all times) as well as the inlet velocity. Under the effect of the periodic boundary, the inlet velocity can be stably maintained at 0.13 *m*/*s* so that the flow is significantly disturbed before it departures the exit boundary without affecting the downstream flow. This further validates the effectiveness of the damping region.

**Figure 3.** Comparison of flow velocities at inlet and outlet boundary at different times at t = 60.0 s, t = 70.0 s, t = 80.0 s and t = 90 s.

For all these reasons, it is possible to confirm the numerical characteristics of the periodic cyclic boundary of the flow.

#### **4. Model Verification**

In order for SPH to work in numerical sink simulations, characterizing the correct flow into the outflow boundary is key to achieve a continuous steady flow cycle. In this study, the flow was simulated using periodic boundary conditions and a damping zone. The position and velocity of the particles were firstly initialized, simulating the state of the current at the beginning of the movement. At the same time, during the interaction between the current and the porous media, the flow velocity was considered to change with the variation of the porosity. According to the experiments in the literature, porosity values of 0.349 and 0.475 [7] are the most commonly used and hence have been considered in this work to verify the model adopted.

The numerical water tank used for verification is shown in Figure 4 below [7]. The numerical flume has a pore structure placed 2.915 m away from the inlet boundary. The numerical model tank has a length of 6.0 m, an inlet flow rate of 0.13 *m*/*s*, a pore structure long 0.15 m and elevated 0.075 m, a sink channel slope of 0.005, and a viscosity of 1.0 <sup>×</sup> <sup>10</sup>−<sup>6</sup> (Pa·s). Finally, the flow field characteristics were simulated at the porosity of 0.349 and 0.475.

**Figure 4.** Schematic diagram of the numerical water tank.

The flow fields of velocity vectors for porosity equal to 0.349 (Figure 5a) and 0.475 (Figure 5b) are shown below in Figure 5.

**Figure 5.** (**a**) represents the velocity field using porosity = 0.349; (**b**) represents the velocity field using porosity = 0.475. U: Longitudinal velocity (*m*/*s*).

It can be seen from Figure 5 that the flow simulation based on the periodic boundary basically simulates the correct natural movement of the water flow. The internal and external flow fields and back vortex flow fields of the porous media are also visible in Figure 4. The flow field characteristics, the flow velocity U, V, and water depth H were processed dimensionless to verify the flow field velocity distribution characteristics of the porous media at 0.475 and 0.349 and results are shown in Figures 6 and 7, respectively.

**Figure 6.** Comparison of longitudinal and vertical velocity for porosity 0.475. Legend: Red open circle: numerical simulation of longitudinal velocity of smooth particle hydrodynamics (SPH); Green dash dot: numerical simulation of longitudinal velocity presented by Chan et al. [7]; Blue solid line: experimental longitudinal velocity values from Chan et al. [7]; Red solid circle: numerical simulation of vertical velocity of SPH; Green solid triangle: numerical simulation of vertical velocity of Chan et al. [7]; Blue solid square: experimental vertical velocity values from Chan et al. [7]. H: water depth (m); U: longitudinal flow rate (*m*/*s*), V: vertical flow rate (*m*/*s*), U0: Inlet boundary velocity (*m*/*s*). X is the distance between the particle and the left boundary of the pore structure, H is the height of the pore structure, and X/H is a dimensionless treatment of the X axis. The pore structure is placed at the origin of this system.

The numerical simulation of the longitudinal and vertical velocity distributions at different locations around the porous media fit well (±5–10%) with the velocity curves measured for the experimental experiments. This confirms that the SPH method is feasible for the simulation of the flow field around the porous media and the analysis of the longitudinal and vertical velocity at a porosity of 0.475.

According to the above model verification displayed in Figure 7, it can also be seen that the SPH method can simulate the flow field inside and around the porous media at a porosity of 0.349, and results are consistent with the experiment of the flow field inside and outside the porous media. The model error is between 10% and 15%.

Observing Figures 6 and 7, which show the average flow direction and vertical velocity components at different locations, it is also possible to notice specific features that are worth to be highlighted. At X/H = 0.33, the velocity gradient is large. In the vertical position, the upstream edge of the structure (X/H = 0.33) is 1 < Y/h < 1.5, and the velocity is small. At the downstream edge of the structure (X/H = 1.67), according to its velocity curve, it can be seen that the permeability of porous media is good, and its reverse velocity is small.

**Figure 7.** Comparison of longitudinal and vertical velocity for porosity 0.349 Legend: Red open circle: numerical simulation of longitudinal velocity of SPH; Green dash dot: numerical simulation of longitudinal velocity presented by Chan et al. [7]; Blue solid line: experimental longitudinal velocity values from Chan et al. [7]; Red solid circle: numerical simulation of vertical velocity of SPH; Green solid triangle: numerical simulation of vertical velocity of Chan et al. [7]; Blue solid square: experimental vertical velocity values from Chan et al. [7]. H: water depth (m); U: longitudinal flow rate (*m*/*s*), V: vertical flow rate (*m*/*s*), U0: Inlet boundary velocity (*m*/*s*). X is the distance between the particle and the left boundary of the pore structure, H is the height of the pore structure, and X/H is a dimensionless treatment of the X axis. The pore structure is placed at the origin of this system.

X/H = 0.33 indicates that the position is 0.02475 m from the left boundary of the pore structure which is placed at the origin. This location is characterized by the velocity distribution of the flow field when the water flow is coupled with the left boundary of the pore structure. Because this is a periodic boundary and a damping region is added to the flow field, the damping region typically causes some disturbance to the flow field in the initial stage, but the effect of the damping region on the flow field velocity gradually decrease after the particles run for 20 s. It was verified within the model that the deviation has a small effect on the flow field around the pores, and it can be seen that the flow field around the pores is not greatly affected by the damping area.

However, in order to ensure that in the oncoming flow field of the pore structure, the flow simulation of the three parts of the pore flow field and the back vortex field of the pore structure have good convergence, the convergence analysis (Figure 8) was completed.

**Figure 8.** Solid purple triangle D = 0.0104; Dashed line when D = 0.01. X/H = −1.4 indicates the oncoming flow field of the pore structure; X/H = 1.67 indicates that the water flows through the pore structure field; X/H = 4.47 indicates the back vortex field of the pore; Y (m) indicates that the particles of the water flow from the horizontal plane vertical distance; U (*m*/*s*) represents the longitudinal water flow velocity; V (*m*/*s*) represents the vertical water flow velocity.

Figure 8 depicts the convergence analysis of the vortex area behind the pore structure and the vortex structure in the oncoming flow area when the porosity is 0.475 and 0.349, respectively. Analyzing Figure 8, it is possible to state that the particles pass through these three areas, and that due to the good convergence and good curve fitting that can be observed, the simulation with different particle sizes still meets the requirements of good convergence.

In practical engineering, the seepage problem of permeable structure will be helpful to the analysis of important engineering structures. It is necessary to study the flow state of fluid in permeable structure. The traditional velocity measurement method is difficult to be used in permeable buildings. Particle image velocimetry (PIV) and particle tracking velocimetry (PTV) are other robust experimental techniques to solve this problem [28–35]. In this study, the Lagrangian particle tracking method is used to measure the velocity distribution along the water depth at different positions.

#### **5. Model Application**

According to the model verification previously described in Section 4, it can be confirmed that the numerical simulation method based on SPH simulates appropriately the water flowing on the porous media, thereby it was decided to apply this model to a water tank and a porous media featuring periodic boundaries and a damping zone. The length of the domain is 60 m, at a porosity value is 0.5, particle size is 0.15, pore length, width, and height ranges within 0.6 m and 1.5 m, water depth considered corresponds to 2.4 m, viscosity is 1.0 <sup>×</sup> 10−<sup>6</sup> (Pa·s), structure coordinates are 30.0, 0.0. The schematic diagram of this configuration is displayed in Figure 9.

**Figure 9.** Schematic diagram of numerical water tank.

#### *5.1. Flow Field Velocity Distribution Diagram of Di*ff*erent Volumes of Porous Media*

Figure 10 shows the flow field velocity distribution of different volumes of porous media at the porosity of 0.5. According to the flow field displayed below, as the volume of the porous media increases (from top to bottom), the longitudinal flow velocity U of the upper layer of the porous media gradually increases, and the upward-incidence flow region is also expanding, but the enlarged region is not well defined. During the expansion process, the back vortex area expands, and multiple vortex motions are formed. The vortex area increases with the volume of the porous media, and begins to form when the vertical section of the porous media is equal to 1.2 m × 1.2 m. While there is an increase of pore structure volume, the back vortex area expands as well but the eddy strength decreases.

**Figure 10.** Flow field velocity distribution of different volumes of porous media at a porosity of 0.5. U: Longitudinal velocity (*m*/*s*). (**a**) Pore structure: 0.6 <sup>×</sup> 0.6 (m2); (**b**) Pore structure: 0.9 <sup>×</sup> 0.9 (m2); (**c**) Pore structure: 1.2 <sup>×</sup> 1.2 (m2); (**d**) Pore structure: 1.5 <sup>×</sup> 1.5 (m2)

#### *5.2. Analysis of Flow Field Inside and Outside the Porous Media*

Figure 11 displays the horizontal and vertical velocity of the flow in the upstream area of 0.6 m × 0.6 m pore structure (the area where water flows through the pore structure and the vortex area). In this Figure 11, as well as for Figures 12–14, x = 25–29.0 m represents the horizontal coordinate range of the upstream area, x = 30.1–30.5 m represents the horizontal coordinate range of the flow through the pore structure, and x = 31.0–36.0 m represents the horizontal coordinate range of the vortex area.

**Figure 11.** Characteristics of longitudinal and vertical flow fields when the vertical section of the pore is 0.6 m × 0.6 m; solid line: x = 25.0; x = 30.1; x = 31.0; dashed line: x = 26.0; x = 30.2; x = 32.0; dotted line: x = 27.0; x = 30.3; x = 33.0. Long dashed line: x = 28.0; x = 30.4; x = 34.0; double dotted line: x = 29.0; x = 30.5; x = 35.0; dotted line: x = 36.0.

Moving from x = 0.0 m towards the pore structure (located at 30.0 m) when the distance between the water particles and the pore material becomes smaller, the longitudinal velocity *u* on the upper side of the pore material increases, as well as the internal velocity of the pore material. The average velocity of this area tends to be stable at about 1 *m*/*s*, and the vertical average velocity V ranges between −0.2–0.1 *m*/*s*.

When x = 30.1–30.5 m, the water flows through the pore area. Increasing x, at the initial part or the pore structure, the longitudinal velocity *u* increases from 0.15 *m*/*s* to 1.5 *m*/*s*, showing the trend of short-term reflux. Therefore, it can be considered that the pore structure of 0.6 m × 0.6 m has an obvious energy dissipation effect, while in the upper water area of the pore structure, the flow velocity shows a slow upward trend and tends to be stable at 1.3 *m*/*s*. In the meantime, the vertical average velocity *V* ranges between −0.2 and 0.4 *m*/*s*.

When x = 31.0–36.0 m, this stage is the eddy current area at the back of the pore structure. It can be seen from the change of the back vortex flow field that with the increase of the volume of porous media, there is a turning point at the speed of H (m) = 0.6 m. There are two stages of velocity curves, before and after the turning point, that can be described independently. The first stage represents the velocity curve inside the pore structure, and the second stage corresponds to the velocity curve outside the pore structure. When *H* is less than the height of pore structure, the velocity of the flow particles increases with the increase of distance from the right boundary of pore structure, with an average velocity of 0.5 *m*/*s*; when the water level *H* is greater than the height of the pore structure, the velocity decreases getting far away from the right boundary of the pore structure, while when x > 30.6 m, and H < 0.6, the velocity of the flow particles on the right boundary of pore structure decreases with the increase of distance from the right boundary of the pore structure. When x = 36.0 m, the flow particle velocity reaches the peak of 1.3 *m*/*s*. For this scenario, the vertical average velocity V is between −0.2 and 0.1 *m*/*s*.

Figure 12 displays the horizontal and vertical velocity of the flow in the upstream area of 0.9 m × 0.9 m pore structure (the area where water flows through the pore structure and the vortex area).

**Figure 12.** Characteristics of longitudinal and vertical flow fields when the vertical section of the pore is 0.9 m × 0.9 m; solid line: x = 25.0; x = 30.1; x = 31.0; dashed line: x = 26.0; x = 30.2; x = 32.0; dotted line: x = 27.0; x = 30.3; x = 33.0. Long dashed line: x = 28.0; x = 30.4; x = 34.0; double dotted line: x = 29.0; x = 30.5; x = 35.0; dotted line: x = 36.0.

When x = 25–29.0 m, with the increase of x, that is, when the distance between water particles and pore structure becomes smaller, the particle flow velocity in the upstream area of pore structure slightly decreases when the longitudinal velocity u on the upper side of pore structure is 0.6 m higher than the one recorded on the side, while the internal velocity growth of the pore structure is still almost negligible. The pore structure is more stable when the overall average velocity is 0.6 m higher than the side length, and the average up-flow velocity tends to be stable when it is close to 1 *m*/*s*, while the vertical average velocity V ranges between −0.2 and 0.1 *m*/*s*. With the increase of x, the longitudinal velocity u on the upper side of the porous material increases slightly, and the decrease observed on the left side is equal to the height of the porous material. The total velocity of the upwelling increases gradually, and tends to be stable near 0.9 *m*/*s*. The vertical velocity of the upwelling area remains almost stable at 0 *m*/*s*.

When x = 30.1–30.5 m, water flows through the pore structure area. With the increase of x, the longitudinal velocity *u* of the fluid in the pore structure decreases gradually, but this phenomenon is not significant, and the reflux trend weakens. Therefore, it can be judged that the energy dissipation effect of the pore structure is slightly stronger when the side length is 0.9 m than when the side length is 0.6 m, while in the upper water area of the pore structure, the flow velocity shows a slow upward trend. When the height h is constant, the flow velocity decreases with the increase of x, and the overall flow velocity still increases with the increase of H.

When x = 31.0–36.0 m, an eddy current area at the back of pore structure can be noticed. It can be seen from the change of the back vortex flow field that there is a turning point in the velocity distribution curve at H = 1.5 m. When H is less than the height of porous material, the velocity increases with the increase of x, and the average velocity is 0.25 *m*/*s*; when h is greater than the height of porous material, the velocity decreases with the increase of x, and reaches 1.7 *m*/*s* at x = 36.0 m, and the back vortex strength decreases with the increase of the distance from the hole. The vertical velocity is relatively stable, in the range of −0.15 and 0.1 *m*/*s*.

Figure 13 displays the horizontal and vertical velocity of the flow in the upstream area of 1.2 m × 1.2 m pore structure (the area where water flows through the pore structure and the vortex area). When x = 25–29.0 m, with the increase of x, the particle velocity in the upstream area of the pore structure decreases slightly when the longitudinal velocity u on the upper side of the pore structure is 0.9 m longer than that on the side, and the average velocity of upwelling tends to be stable at the place close to 0.8 *m*/*s*, and the vertical average velocity V is between −0.15–0.1 *m*/*s*. When x = 30.1–30.5 m, water flows through the pore area.

**Figure 13.** Characteristics of longitudinal and vertical flow fields when the vertical section of the pore is 1.2 m × 1.2 m; solid line: x = 25.0; x = 30.1; x = 31.0; dashed line: x = 26.0; x = 30.2; x = 32.0; dotted line: x = 27.0; x = 30.3; x = 33.0. Long dashed line: x = 28.0; x = 30.4; x = 34.0; double dotted line: x = 29.0; x = 30.5; x = 35.0; dotted line: x = 36.0.

With the increase of x, the longitudinal velocity U of the fluid in the pore structure is gradually stable at 0.1 *m*/*s*, and the reflux trend is not strong. It can be seen that for 1.2 m × 1.2 m porous material, the energy dissipation effect of porous material is also stronger than before. When x = 31.0–36.0 m, when H is greater than the hole height, the velocity decreases with the increase of x, and the velocity at the right boundary of the pore structure decreases with the increase of the distance from the right boundary of the hole to 2.0 *m*/*s*. The back vortex intensity also decreases with the increase of x. The vertical velocity is relatively stable, in the range of −0.20 and 0.35 *m*/*s*.

Finally, Figure 14 displays the horizontal and vertical velocity of the flow in the upstream area of 1.5 m × 1.5 m pore structure (the area where water flows through the pore structure and the vortex area). When x = 25–29.0 m, the overall velocity of the upwelling increases gradually, and tends to be stable (approximately 0.75 *m*/*s*), while the vertical velocity ranges between −0.2 and 0.2 *m*/*s*.

When x = 30.1–30.5 m, water flows through the pore area. With the increase of x, the longitudinal velocity *u* of the fluid in the pore structure tends to 0, and the reflux trend is not noticeable. It can be concluded that the hole has obvious energy dissipation effect on the hole of 1.5 m × 1.5 m, while in the upper water area of the hole, the velocity keeps rising slowly. With the increase of distance, the water content of the upper part of the pore structure shows a slowly rising trend, the vertical velocity increases at the same time. When x = 31.0–36.0 m, when h is less than the height of porous material, the flow velocity increases with the increase of the right boundary distance of porous material, and the average flow velocity is approximately 0.02 *m*/*s*; when h is greater than the height of porous material, the flow velocity decreases with the increase of the right boundary distance of porous material, and the flow velocity at the right boundary of porous material reaches 2.2 *m*/*s*, and the back vortex strength decreases with the increase of the distance from the hole. The vertical velocity is relatively stable, in the range between −0.35 and 0.35 *m*/*s*.

**Figure 14.** Characteristics of longitudinal and vertical flow fields when the longitudinal section of the pore is 1.5 m × 1.5 m; solid line: x = 25.0; x = 30.1; x = 31.0; dashed line: x = 26.0; x = 30.2; x = 32.0; dotted line: x = 27.0; x = 30.3; x = 33.0. Long dashed line: x = 28.0; x = 30.4; x = 34.0; double dotted line: x = 29.0; x = 30.5; x = 35.0; dotted line: x = 36.0.

In summary, based on the changes of up-flow field tested, it can be stated that by increasing the pore volume, the up-flow velocity typically decreases and the energy dissipation becomes higher, but changes associated with the vertical velocity are not significant with the change of volume for the pore structure. Focusing on the change of the flow field inside and outside the hole it is possible to confirm that the internal velocity of the pore is mainly in the range of 0 and 0.2 *m*/*s*. With the increase of the porous media' volume, the internal velocity difference of the pore decreases gradually, while the particle velocity directly above the pore structure increases with the increase of the pore volume, and the velocity difference of the upper layer flow also gradually increases.

Furthermore, it can be seen from the change of the back vortex flow field that, with the increase of the volume of porous media, there are inflexion points at H (m) = 0.6, 0.9, 1.2 and 1.5 for the four groups of horizontal velocity, respectively, thus forming multiple velocity curves before and after the inflexion point. When H is less than the height of the hole, the flow velocity increases with the increase of x, however when H is greater than the height of the hole, the flow velocity increases with the increase of x. The intensity of the back vortex also decreases with the increase of x. From the change of vertical velocity, it can be observed that the velocity basically varies from −0.5 *m*/*s* to 0.5 *m*/*s*.

#### *5.3. Longitudinal and Vertical Flow Field Distribution under Di*ff*erent Porosity*

Porosity is also an important factor influencing the characteristics of the flow field in the porous media [36]. Therefore, the pore size was 1.2 × 1.2 (m), and the oncoming flow and reflux vortices were analyzed for porosity values of 0.5, 0.25 and 0.125, respectively. The velocity distribution and the flow field along the water depth were analyzed for each porosity investigated.

Figure 15 displays the flow field in the upstream area of pore structure under different porosities 0.5, 0.25, and 0.125 for x = 25.0–29.0. As the porosity decreases, it can be seen that the longitudinal velocity curve of the upstream flow tends to be gradually stable.

**Figure 15.** Longitudinal velocity distribution of the flow field in the upstream area of the pore structure when the porosity is 0.5, 0.25 and 0.125 and x = 25.0–29.0 m. (**a**) nw = 0.5; (**b**) nw = 0.25; (**c**) nw = 0.125.

Figure 16 shows the flow field when the water flows through the pore structure under different porosities 0.5, 0.25 and 0.125 for x = 30.1–30.5. The longitudinal velocity difference of the internal flow field of the pore structure increases with the decrease of porosity.

**Figure 16.** Longitudinal velocity distribution of the flow field inside and directly above the pore structure when the porosity is 0.5, 0.25 and 0.125 and x = 30.1–30.5.0 m. (**a**) nw = 0.5; (**b**) nw = 0.25; (**c**) nw = 0.125.

Figure 17 presents the flow field in the back vortex region with different porosities 0.5, 0.25 and 0.125 for x = 31.0–36.0. With the decrease of porosity, the flow rate changes are not significant, but the back vortex region increases slightly, and the flow rate at the inflection point also initially increases but later decreases in combination with the increase of porosity. When the porosity is 0.25, the particle longitudinal velocity at the inflection point is the maximum calculated.

**Figure 17.** Longitudinal velocity distribution in the back vortex flow field when the porosity is 0.5, 0.25, and 0.125 and x = 31.0–36.0 m. (**a**) nw = 0.5; (**b**) nw = 0.25; (**c**) nw = 0.125.

Figure 18 displays the vertical velocity distribution of the flow field in the upstream area of the pore structure under different porosities 0.5, 0.25, and 0.125 for x = 25.0–29.0. With the gradual decrease of porosity, it can be seen that the vertical velocity difference of each position of the upstream flow is gradually reduced.

**Figure 18.** Vertical velocity distribution of the flow field in the upstream area of the pore structure when the porosity is 0.5, 0.25 and 0.125 and x = 25.0–29.0 m. (**a**) nw = 0.5; (**b**) nw = 0.25; (**c**) nw = 0.125.

Figure 19 shows the vertical velocity distribution of the flow field inside and directly above the pore structure under different porosities 0.5, 0.25, and 0.125 for x = 30.1–30.5. With the gradual decrease of porosity, it can be seen that in the fluid solid interface area, the change of horizontal velocity becomes more evident, and the velocity reduces.

**Figure 19.** Vertical velocity distribution of the flow field inside and directly above the pore structure when the porosity is 0.5, 0.25 and 0.125 and x = 30.1–30.5 m. (**a**) nw = 0.5; (**b**) nw = 0.25; (**c**) nw = 0.125.

Figure 20 presents the vertical velocity distribution of the back vortex flow field under different porosities 0.5, 0.25 and 0.125 for x = 31.0–36.0. As the porosity decreases gradually, it can be seen that when the height H is constant, the vertical velocity difference at each position is characterized by a gradually increasing trend.

**Figure 20.** Vertical velocity distribution in the back vortex flow field when the porosity is 0.5, 0.25 and 0.12 and x = 31.0–36.0 m. (**a**) nw = 0.5; (**b**) nw = 0.25; (**c**) nw = 0.125.

In summary, it can be seen from Figures 15–20 that as the porosity increases, the longitudinal flow velocity of the oncoming flow slightly changes, but the overall change is not significant, plus the longitudinal flow velocity of the back vortex decreases. Considering the vertical flow velocity of the upstream area, despite the increase of porosity, this parameter remains within the range between −0.2 *m*/*s* and 0.2 *m*/*s*, and the vertical flow velocity through the porous media decreases with the increase of porosity. From the vertical flow velocity of the back vortex, as the porosity increases, the vertical flow velocity tends to decrease slightly overall.

Due to the infiltration of the water, the permeable structure has a significant impact on its discharge capacity [37]. It can be seen from the horizontal and vertical comparison of the flow field in Figure 16 that with the increase of porosity, the flow resistance decreases. The distribution of longitudinal velocity of porous material changes in the upper layer of porous material.

The maximum turbulence intensity near the surface and behind the structure increases proportional to the porosity [38] at the same downstream position.

#### *5.4. Convergence Verification of Pore Logistics Field Model*

In order to verify the convergence of the SPH in the numerical model, the flow field characteristics of each section previously investigated were analyzed by applying multiple particle size distributions [8]. Up-flow area, pore flow area, back vortex area, and the horizontal and vertical velocity distribution along the depth were calculated with the particle size D of 0.096 m, 0.01 m, and 0.104 m as shown in Figures 21 and 22 to verify convergence, at the selected positions x = 28.0 m, 30.3 m, 33 m.

**Figure 21.** Longitudinal velocity distribution of the inflow field with different particle spacing (**a**), longitudinal velocity distribution of the sea current flowing through the porous media with different particle spacing (**b**), and longitudinal velocity distribution of the rear vortex field with different particle spacing (**c**).

**Figure 22.** Vertical velocity distribution of the inflow field with different particle spacing (**a**), vertical velocity distribution of the sea current flowing through the porous media with different particle spacing (**b**), and vertical velocity distribution of the rear vortex field with different particle spacing (**c**).

It was previously noticed that when the height of the pores is 1.2 m, the flow direction is opposite to the others when D = 0.096. This elevation corresponds to the upper flow area of the pore structure and justify the presence of turbulence as shown in the graph. However, results displayed in Figures 21 and 22 show that the velocity distributions of the internal and external flow fields formed by the oncoming flow field, the back vortex flow field and the water flow through the porous media are consistent despite different particle spacing and the error range is between 10% and 15%, hence providing an acceptable grade of convergence.

Finally, Figures 23 and 24 show the back vortex field of the pore structure when the porosity is 0.349 and 0.475, respectively.

**Figure 23.** Back vortex field of the pore structure when the porosity is 0.349. U: Longitudinal velocity (*m*/*s*).

**Figure 24.** Back vortex field of the pore structure when the porosity is 0.475. U: Longitudinal velocity (*m*/*s*).

It can be seen from the flow field distribution diagrams in Figures 23 and 24 that the energy in the vicinity of the pore structure presents a dissipative trend and shows different vorticities under different porosities. When the porosity is 0.475, the energy dissipation in the vicinity of the pore structure is weak, the vortex flow pattern is more noticeable and multiple vortex flow patterns are formed on the back of the pore structure. However, when the porosity is 0.349, the energy dissipation near the pore structure is relatively strong, and the intensity of the vortex flow is relatively weak.

To enable academic colleagues to replicate the presented numerical data and have a better estimation of the numerical capabilities of the numerical method, details about the total number of particles used, the particle size, the physical simulation times plus the CPU time and the CPU cores used are displayed and summarized in Table 1 for all the configurations tested (D = 0.0104 m, D = 0.01 m, D = 0.096 m, D = 0.1 m, and D = 0.104 m).


**Table 1.** Simulation details for each configuration tested.

#### **6. Conclusions**

Nowadays, a variety of man-made structures such as submerged breakwaters, fishing reefs, outfall protections, armor layers, rubble-mound or berm breakwaters, built of gravel or artificial units, can be found in coastal areas worldwide for multiple purposes (e.g., coastal flooding protection, coastal erosion protection). All these structures are characterized by the fact that either some of their layers or their full structure are porous, hence they are permeable to flows induced by waves and currents. Therefore, the hydrodynamics, the stability and performance of these structures are dependent on the characteristics of the waves and their interaction with permeable material and it is crucial to provide mathematical formulations to model these features and the role played by the porous material considering its geometry, its location, and its characteristics (which could influence the wave propagation, diffusion, overtopping, and the consequent turbulence generated).

This paper enabled the calibration and validation of a SPH model utilized to characterize the influence of porous media on flow propagation and outcomes can be summarized as follows:

1. Based on the analysis of the area upstream the porous media, the longitudinal flow velocity of the ascending flow decreases slightly with the increase of the volume of the pore structure;

2. Based on the analysis of the internal and external flow fields of the porous media, the internal flow velocity of the porous media for the configurations tested is in the range of 0–0.2 *m*/*s*. It was possible to notice that with the increase of the volume of the porous media, the flow velocity of the upper layer of the porous media firstly tends to accelerate and then stabilizes. The longitudinal flow velocity of the upper layer of the flow increases with the increase of water depth, while the vertical velocity fluctuates sharply at the fluid-solid boundary;

3. Based on the back vortex analysis, there is an inflection point in the flow velocity distribution of the upper and lower layers of the porous media. As the distance from the porous media increases, the flow velocity increases, so that the velocity distribution curve of the back vortex flow forms different concavities and convexities before and after the porous media, and as the pore volume gradually increases, a plurality of vortices are formed on the backflow surface. With the rotary motion, the area of the back vortex increases as the volume of the pore increases, but the strength decreases as the volume increases.

4. Overall, with the increase of porosity, the longitudinal flow and back vortex flow slightly decrease in the vertical direction, and the SPH-based model calculation has good convergence (error between experimental and numerical results within 10%).

**Author Contributions:** Conceptualization, Q.G.; methodology, S.W. and Q.G.; software, S.W. and Q.G.; validation, S.W. and Q.G. and M.R.; formal analysis, S.W. and Q.G.; investigation, S.W., Q.G. and M.R.; resources, S.W. and Q.G.; data curation, S.W., Q.G. and M.R.; writing—original draft preparation, S.W., Q.G. and M.R.; writing—review and editing, S.W., Q.G. and M.R.; visualization, S.W. and Q.G.; supervision, Q.G. and M.R.; project administration, Q.G.; funding acquisition, Q.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Zhejiang Province Natural Science Foundation, grant number LQ16E090001, the National Natural Science Foundation of China, grant number 51509134.

**Acknowledgments:** The author is very grateful to Ningbo University, the National Natural Science Foundation of China and the Natural Science Foundation of Zhejiang Province. All authors are very grateful to Dr Songdong Shao for his constructive comments and continuous inspiration to complete this work.

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

#### **References**


© 2020 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* **An Experimental Study of Focusing Wave Generation with Improved Wave Amplitude Spectra**

**Guochun Xu 1,\*, Hongbin Hao 2, Qingwei Ma <sup>2</sup> and Qinqin Gui <sup>1</sup>**


Received: 15 October 2019; Accepted: 26 November 2019; Published: 28 November 2019

**Abstract:** We experimentally investigate the generating results of space-time focusing waves based on two new wave spectra, i.e., the quasi constant wave amplitude spectrum (QCWA) and the quasi constant wave steepness spectrum (QCWS), in which amplitude and steepness for each wave component can be adjusted with fixed wave energy. The wavemaker signal consists of a theoretical wavemaker motion signal and two different auxiliary functions at two ends of the signal. By testing a series of focusing waves in a physical wave tank, we found that with given wave energy, the QCWA spectrum can produce a focusing wave with larger crest elevation and farther focusing location from the wavemaker flap, as compared with the QCWS spectrum. However, both spectra lead to larger focusing wave crests when the wave frequency bandwidth was narrowed down and a positive correlation between the generated relative wave crest elevation and the input wave elevation parameter. The two spectra produce different focusing wave positions for the same wave frequency range. We also found that the focusing time strongly relates to the energy of the highest-frequency wave component of the wave spectrum.

**Keywords:** focusing waves; wave amplitude spectra; space-time parameter; experimental investigations

#### **1. Introduction**

Focusing wave is a special water wave, different from regular wave or stochastic wave, with a single large wave crest when it happens [1,2]. On the basis of investigations of the triggering mechanisms of focusing waves [3–6], many causes have been identified, such as the space-time focusing of transient waves, wave–current interaction, geometrical focusing due to the seabed topography, atmospheric forcing, nonlinear self-focusing, etc. If the wave heights of focusing waves exceed two to 2.2 times their significant wave heights, they are generally defined as freak waves or rogue waves [7,8]. Hence, focusing waves in laboratory are often employed to model freak wave events observed in extreme sea state [9,10], in order to better understand the generation process, the mechanisms of those extreme waves, and the hydrodynamic loads on floating or fixed ocean structures in extreme sea environments [11–13].

Over the past few decades, various theoretical models, such as wave energy focusing, wave–wave interaction, wave-air coupling, etc. [14–16], have been established to express the generation mechanisms of focusing waves. Among them, the space-time focusing theory of wave energy and wave modulation instability (also named nonlinear self-focusing) theory are the most widely used and are often applied to numerical simulations of focusing or extreme waves. The former assumes that the focusing wave consists of a number of small harmonic wave components which can be superposed to form the focusing wave. Research based on the space-time focusing theory was initially carried out by the linear wave superposition [17,18] and, then, it was followed by low-order wave–wave interactions [19], and wave directional spread [20]. Recently, based on spatio-temporal focusing of wave energy, focusing waves have been produced in fully nonlinear potential or viscous numerical wave tanks [21–24]. As for the latter, i.e., wave modulation instability, a wave group is designed to be composed of carrier waves and their sideband waves. When the sideband disturbance occurs in a wave travelling process, the wave energy of carrier waves is transferred to their sideband waves, leading to wave energy focusing. Waves generated this way are are also called rogue waves or freak waves. Theoretical study on wave modulation instability dates back to the last sixties with investigations on a Stokes wave train with small perturbation [25]. Then, studies were extended to nonlinear four-wave interactions and random wave group [26,27]. In the meantime, several mathematical equations were derived under different physical assumptions, such as the nonlinear Schrödinger equation, the Davey–Stewartson system, the Korteweg–de Vries equation, and the Kadomtsev–Petviashvili equation. Detailed reviews of these mathematical models can be found in [3,28].

On the basis of the abovementioned generating mechanisms and theoretical models, the focusing waves have been mainly produced by two methods in physical wave flumes. One is to produce focusing waves through space-time focusing of wave energy. To achieve wave energy focusing at a specified position and time, the phases of wave components are modulated as zero or π/2. This phase modulation method has been adopted in the single wave model [19], the double wave model [29], the triple wave model [30], and the NewWave model [31]. The main difference among those models is the ways to specify the amplitudes of wave components. In the first three models, their amplitudes are assigned by some predefined wave spectrum, such as the constant wave amplitude spectrum (CWA), the constant wave steepness spectrum (CWS), and the random sea wave spectrum or their combined spectra, whereas the wave component amplitudes in the NewWave model are determined by the autocorrelation function of the wave energy density spectrum. In addition to the above phase modulation method, focusing waves can also be spatio-temporally produced by means of the wave dispersion method, i.e., focusing wave components are individually generated and their frequencies linearly vary with the largest one at the start [32]. By contrast, the phase modulation method can produce focusing waves with higher frequency wave components. When the wave dispersion method is adopted to generate focusing waves, the higher frequency waves are severely constrained by the stroke limitation of wavemaker.

The other way to generate focusing waves, in a laboratory, is to employ wave modulation instability. On the basis of this mechanism, Li et al. [33] experimentally observed focusing wave occurrence in a random wave group by adjusting the wave steepness parameter and the Benjamin–Feir Index (BFI). By assigning carrier wave amplitude and steepness, Chabchoub et al. [34] generated super focusing waves in their physical wave flume. Moreover, several breather solutions of the nonlinear Schrödinger equation, such as Kuznetsov–Ma solution, Akhmediev breather, and the Peregrine breather solution, are sometimes adopted to experimentally simulate extreme waves in a physical wave tank [35,36]. Apart from above methods, recently a phase-amplitude iteration scheme based on space-time focusing of wave energy has been developed to generate tailored focusing wave [37,38]. Nevertheless, the investigation from Deng et al. [39] indicated that the above phase-amplitude iteration scheme may be less effective for focusing waves with severe phase coupling. The phase-amplitude iteration method probably can be improved by directly adjusting wave component amplitude or steepness.

In order to better understand focusing waves generated under various experimental conditions, several experimental investigations have also been carried out. The testing results from Liu et al. [40] showed that some additional high-order wave components are produced, and the amplitudes of these extra components increase when focusing wave elevation becomes large. The frequency parameter of focusing wave was experimentally examined by Ma et al. [41] who demonstrated that focusing waves with a wider frequency range transfer more wave energy to their high-frequency components. The local wave steepness of focusing waves with three wave amplitude spectra types, (i.e., the linear wave steepness spectrum, the CWA spectrum, and the CWS spectrum) were examined by Wu and Yao [42] in a physical wave-current flume. Their results demonstrate that focusing waves generated by the linear steepness wave spectrum have larger local wave steepness. In addition, focusing waves

with the CWA spectrum were investigated by Baldock et al. [19] who observed that the focusing position, time, and wave crest elevation increase along with a larger amplitude parameter. However, despite all the wave amplitude spectra used to generate focusing waves in published researches, their wave component amplitudes or steepness are not adjustable under assigned wave energy. Indeed, in the recently developed phase-amplitude iteration scheme [37] used for focusing wave generation, the steepness of wave components are fixed once the wave component amplitudes are determined. This means that the wave nonlinearity effect arising from phase coupling among wave components cannot be better controlled in focusing wave generation. Thus, more flexible spectra with adjustable amplitude and steepness of the components should be explored in the laboratory, in order to better represent real focusing waves or freak waves with strong nonlinearity [43].

In this study, two new wave spectra, i.e., quasi constant wave amplitude spectrum (QCWA) and quasi constant wave steepness spectrum (QCWS), are developed by modifying the conventional CWA and CWS spectra [44] used to produce focusing waves. In comparison to previous wave spectra, the advantage of the new wave spectra is that the amplitudes and steepness of wave components can be adjusted by adjusting the water depth of wave tank. On the basis of the two new wave spectra, two-dimensional nonbreaking focusing waves are spatio-temporally generated in a physical wave tank and the parameters of the generated waves are investigated under different component groups.

#### **2. Generating Principle of Focusing Waves in a Physical Wave Tank**

#### *2.1. Experimental Setup*

The experiment on focusing wave generation was carried out in the towing wave tank at Harbin Engineering University. The experimental sketch is shown in Figure 1. The wave tank has a length of 108 m, a width of 7 m, and a depth of 3.5 m. The hydraulic flap-type wavemaker is installed at one side of the wave tank and generates regular waves and random waves, with the maximum wave height of 0.4 m and the wave period ranging from 0.4 s to 4.0 s. At the other side, an absorbing shore is arranged to reduce the reflective waves from the end wall of the wave tank.

**Figure 1.** Experimental arrangement sketch.

To record the free surface elevations, 24 wave gauges were installed with an interval of 0.4 m along the length of the wave tank. When the space parameter of focusing waves was set as 30 m, 40 m, and 50 m, the gauge array could be shifted with the first wave gauge (the far left one in Figure 1) installed at 21.905 m, 31.495 m, and 38.905 m from the flap, respectively. The hydraulic piston movement was also monitored by optic equipment (including markers and the matching cameras) to validate the accuracy of output signals.

The experimental procedure is summarized in Figure 2. To generate the desired focusing waves in the physical wave tank, it is essential to have appropriate and practical wave-making signals. The theory and method used in this research for signal generation are described in the following Sections 2.2 and 2.3.

**Figure 2.** Flow chart of focusing wave generation.

#### *2.2. Theoretical Wave-Making Signal*

Focusing waves are generated by the spatio-temporal focusing principle of wave energy. The origin of the coordinate system is defined at the intersection point of the flap and the still water surface. The z-axis is vertically upward, and the x-axis is horizontally towards the wave travelling direction. The free surface of the focusing wave in the defined coordinate system can be expressed by:

$$\eta = \sum\_{n=1}^{N} \mathbf{a}\_{\rm n} \cos[\mathbf{k}\_{\rm n}(\mathbf{x} - \mathbf{x}\_{\rm f}) - \omega\_{\rm n}(\mathbf{t} - \mathbf{t}\_{\rm f})] \tag{1}$$

where xf and tf are focusing location and time and ω<sup>n</sup> is the frequency of each wave component, which linearly increases from ω<sup>1</sup> to ωN, as expressed by Equation (2). In this work, an is newly developed to be a variable which is determined by wave spectra, i.e., the quasi constant wave amplitude spectrum and the quasi constant wave steepness spectrum. They can be expressed by Equations (3) and (4), respectively.

$$
\omega\_{\mathbf{n}} = \omega\_1 + (\mathbf{n} - 1) \times \frac{\omega\_{\mathbf{N}} - \omega\_1}{\mathbf{N} - 1} \tag{2}
$$

$$\mathbf{a}\_{\mathbf{n}} = \frac{\mathbf{A}\_{\mathbf{f}}}{\mathbf{N}} \frac{\cosh(\mathbf{k}\_{\mathbf{n}} \mathbf{h})}{\sinh(\mathbf{k}\_{\mathbf{n}} \mathbf{h})} \tag{3}$$

$$\mathbf{a}\_{\mathbf{n}} = \frac{\mathbf{A}\_{\mathbf{f}}}{\mathbf{k}\_{\mathbf{n}} \sum\_{\mathbf{n}=1}^{N} 1/\mathbf{k}\_{\mathbf{n}}} \frac{\cosh(\mathbf{k}\_{\mathbf{n}} \mathbf{h})}{\sinh(\mathbf{k}\_{\mathbf{n}} \mathbf{h})} \tag{4}$$

From Equations (3) and (4), it can be seen that cosh(knh)/sinh(knh) approaches one when the water depth tends to be infinite. Thus, the previous CWA and CWS spectra used to generate focusing waves are two special expressions of the QCWA and QCWS spectra at infinite water depth.

The wave component amplitudes of the QCWA spectrum and wave steepness of the QCWS spectrum at different water depths are demonstrated in Figure 3a,b, respectively. Figure 3a also shows that the wave component amplitude of the QCWA spectrum approaches that of the CWA spectrum when water depth increases. For a certain water depth, the wave component amplitude gradually becomes a constant as the wave frequency becomes higher. This is because these high-frequency wave components have short wave lengths and can be considered as deep-water waves. Therefore, Equation (3) expresses the CWA spectrum for high frequency range, i.e., an = Af/N. The similar correlation could also be seen between the QCWS spectrum and the CWS spectrum in Figure 3b. According to QCWA (or QCWS) spectrum, for the finite water depth, the crest elevation of the generated focusing wave (defined as Af 0) is equal to an (calculated by the linear wave theory), which is larger than Af resulting from the CWA (or CWS) spectrum. According to the focusing wave-free surface elevation expressed by Equation (1), the hydraulic piston movement to drive wavemaker flap yields

$$\mathbf{s}(\mathbf{t}) = \sum\_{\mathbf{n}=1}^{N} \frac{\mathbf{a}\_{\mathbf{n}}}{\mathbf{F}\_{\mathbf{n}}} \cos[\omega\_{\mathbf{n}} \mathbf{t} + (\mathbf{k}\_{\mathbf{n}} \mathbf{x}\_{\mathbf{f}} - \omega\_{\mathbf{n}} \mathbf{t}\_{\mathbf{f}})] \tag{5}$$

where,

$$\mathbf{F}\_{\mathbf{n}} = \frac{4\mu\_{\mathbf{n}}^2}{\mathbf{g}\mathbf{H}\_1} \frac{\cosh(\mathbf{k}\_\mathbf{n}\mathbf{h}) \left[\cosh(\mathbf{k}\_\mathbf{n}(\mathbf{h} - \mathbf{H}\_0)) - \cosh(\mathbf{k}\_\mathbf{n}\mathbf{h}) + \mathbf{k}\_\mathbf{n}\mathbf{h}\sinh(\mathbf{k}\_\mathbf{n}\mathbf{h})\right]}{\mathbf{k}\_\mathbf{n}^2(2\mathbf{k}\_\mathbf{n}\mathbf{h} + \sinh(2\mathbf{k}\_\mathbf{n}\mathbf{h}))} \tag{6}$$

and H1 = 2.3 m for the wavemaker of this experiment. The displacement signal calculated by Equation (5) is the theoretical wave-making signal to generate focusing waves.

**Figure 3.** Two wave spectra changing with water depth: (**a**) The quasi constant wave amplitude QCWA spectrum and (**b**) the quasi constant wave steepness QCWS spectrum.

#### *2.3. Implementing Wave-Making Signal*

The theoretical wave-making signal cannot be directly employed to generate focusing waves, because of nonzero displacement at initial and terminating instants. Practically, the wavemaker flap needs to slowly start from its static status, and it has to gradually terminate at static status after the wave generation is completed. Thus, the theoretical wave-making signal needs to be processed to match the actual motion of the flap. An example of the theoretical signal of a QCWA focusing wave is illustrated by the red line in Figure 4 with Af = 0.1 m, xf = 50 m, tf = 25 s, N = 32, and ω<sup>n</sup> = 1.336 rad/s~2.695 rad/s, and the wave component amplitude is computed by Equation (3).

**Figure 4.** Theoretical signal and processed signal by ramp function.

From the theoretical wave-making signal, it is seen that the initial displacement is smaller than zero and the terminating value is larger than zero. This is generally remedied by multiplying the theoretical signal by a ramp function, as shown in Figure 4. It is observed that the flap gradually starts from the static position and stops at the initial static position, matching with the actual motion requirements of wavemaker.

However, by comparing the rump function modified and the original wave-making signals in Figure 4, we observe that the theoretical movement signal has been distorted by the ramp function at the beginning and terminating stages. Considering the strong transient characteristics of focusing waves, the ramp function method may not be suitable for processing the theoretical wave-making signal. In this study, a new scheme is proposed to generate the practicable wave-making signal of focusing waves through inserting a small section auxiliary signal expressed as

$$\mathbf{s}\_{0}(\mathbf{t}) = \left(\sum\_{j=1}^{m} \mu\_{j}(\mathbf{t}) \chi\_{j}\right) \times \left(\sin^{\gamma}\left(\frac{\pi}{2\mathbf{t}\_{0}}\mathbf{t}\right)\right) \left(0 \le \mathbf{t} \le \mathbf{t}\_{0}\right) \tag{7}$$

where, the base function μ<sup>j</sup> is written as [1, t, t2] for m = 3 and the coefficient matrix χ could be solved from Equation (8). For matching with the initially static state of the flap, the parameter γ should be larger than one. In this study, γ = 2, t0 = 1.5 s.

$$\begin{cases} \mathbf{s}\_0(\mathbf{t}\_0) = \mathbf{s}(0) \\ \mathbf{s}'\_0(\mathbf{t}\_0) = \mathbf{s}'(0) \\ \mathbf{s}''\_0(\mathbf{t}\_0) = \mathbf{s}''(0) \end{cases} \tag{8}$$

To make the wavemaker flap return to its initial position at the terminating stage, the theoretical stroke signal is extended, and the extension is processed by a ramp down function. The newly proposed scheme is applied to the above focusing wave example in Figure 4, and the modified wave-making signal result is illustrated in Figure 5.

**Figure 5.** Wave-making signal generated by the new scheme.

It is observed that the wave-making signal processed by this research meets the actually wave-making requirements by extensions at two ends of the theoretical signal. This method will be adopted in the following focusing wave generation.

#### **3. Experimental Results and Discussions**

All the tested focusing wave cases in this research are listed in the Appendix A. In order to examine the accuracy of the wavemaker movement controlling system, the hydraulic piston motion of wavemaker is captured and is compared with the input stroke signal. The motion signal comparison of one focusing wave case (i.e., A101f61X50t36 in Appendix A) is plotted in Figure 6, showing that the tested signal has good agreement with the input signal. The maximum relative error is less than 2%.

**Figure 6.** Comparison of input and tested stroke signals.

The repeatability of the whole experiment system also has been examined by each wave gauge recording result in the wave tank. All focusing wave cases are repeatedly tested twice and the percentage error of the maximum wave crest elevation (or the minimum wave trough) at each tested position is estimated by

$$\mathbf{E}\_{\mathbf{r}}^{+(-)} = \left[ \frac{\mathfrak{n}\_1^{+(-)} - \mathfrak{n}\_2^{+(-)}}{\mathfrak{n}\_1^{+(-)}} \right] \times 100\% \tag{9}$$

The relative percentage errors of one focusing wave case (i.e., A101f61X50t36 in Appendix A) are illustrated in Figure 7a,b, showing the maximum errors of the tested wave crests and wave troughs are about 4% and 3%, respectively. The rest of the experimental focusing wave cases are also validated by repeating tests and their maximums of Er <sup>+</sup> and Er − are both less than 5% [45].

**Figure 7.** Experimental repeatability error: (**a**) Relative error of focusing wave crests and (**b**) relative error of focusing wave troughs.

Moreover, the wave modulation instability of the experimental focusing wave cases in this research has also been examined depending on the Benjamin–Feir Index (BFI). The BFI is evaluated by half-frequency width at half maximum of the wave spectrum [46,47]. The corresponding results of all experimental focusing wave cases are summarized in Table A1. We found that the BFIs of the focusing waves are all less than one, which indicates the wave focusing in this research is dominated by the space-time focusing of wave energy, rather than the modulation instability. In the following sections, the effects of the spectrum types and focusing wave parameters will be investigated according to the tested focusing waves. It should be noted that the focusing wave in this study is defined as the wave with the maximum wave crest elevation.

#### *3.1. Wavefree Surface Evolution for QCWA and QCWS Focusing Waves*

To compare the evolving processes of the focusing waves generated by QCWA and QCWS spectra, a focusing wave case with parameters of Af = 0.1 m, tf = 36 s, xf = 50 m, ω<sup>n</sup> = 0.997 rad/s~3.696 rad/s, and N = 32 was tested. The wave amplitude spectra and wave steepness spectra are compared in Figure 8a,b, respectively. It is seen that, the QCWA spectrum provides larger wave amplitudes and steepness for the high-frequency wave components than the QCWS spectrum. This implies the focusing wave generated by the QCWA spectrum will have stronger nonlinearity.

**Figure 8.** Comparison of two wave spectra: (**a**)Wave amplitude spectrum and (**b**) wave steepness spectrum.

The wave-free surface time histories of focusing waves based on the above two spectra at different positions are plotted in Figure 9a,b showing that the two focusing waves both start from a deep wave trough, followed by the largest wave crest to form the focusing wave. The deep trough then arises and disappears again to complete the wave energy focusing and diffusing [48]. The focusing waves based on the two wave spectra both appear at 45.305 m from the flap, being smaller than the assigned position of 50 m. The shift of the focusing point towards upstream was also observed in the study of the CWA focusing waves [49], while the opposite phenomenon was found in investigations on the CWA and CWS focusing waves by Li and Liu [50]. Hennig and Schmittner [51] pointed out that the shift of focusing point is due to the wave group celerity alteration in focusing wave generation.

**Figure 9.** Evolution of two focusing waves: (**a**) The QCWA focusing wave and (**b**) the QCWS focusing wave.

Relative to the focusing wave generating instant in their time history curves (the wave-free surface curves of the wave-free plotted by red line), the two focusing wave-free surfaces are both asymmetry, which implies the focusing wave generation shown in Figure 9 is not only a linear superposition process. Comparing time history curves at each testing position in Figure 9, we observe that the QCWA focusing wave-free surface is steeper and deeper, whereas the QCWS focusing wave-free surface shape is flatter and shallower. This is mainly introduced by their wave amplitude spectra differences, as shown in Figure 8. The QCWS focusing wave has more wave energy in low frequency range which determines the overall shape of focusing waves, while the QCWA focusing wave distributes more wave energy in high frequencies which affects the local wave-free surface.

To further analyze focusing wave results from the aspect of the spectrum, the wave component amplitude spectra are calculated via a fast Fourier transform (FFT) analysis based on recording results of wave gauges at different positions. Because of the disturbance of reflective waves, only a limited time length of wave surface recording is available in FFT. It should be noted that the resolution of wave spectra calculated by FFT is lower than the input wave spectra. From the FFT results shown in Figures 10 and 11, it is clearly seen that the extra lower-frequency (<0.997 rad/s) and higher-frequency harmonics (>3.696 rad/s) are generated in generating processes of focusing waves. The higher-frequency components are mainly the second-order harmonics (3.696 rad/s~7 rad/s). These higher-frequency harmonics of the QCWA spectrum continually become larger and the wave components in input frequency range reduce as the wave group approaches the focus position. Nevertheless, one opposite result is found when the wave group passes the focus position. For the wave spectra variation of the QCWS focusing wave shown in Figure 11, it is observed that the apparent variation of wave component amplitudes occurs within the input frequency range (i.e., the wave energy is redistributed in input wave frequency).

**Figure 10.** Wave amplitude spectra variation of the QCWA focusing wave at different positions.

**Figure 11.** Wave amplitude spectra variation of the QCWS focusing wave at different positions.

Moreover, for qualitatively identifying the nonlinear effects on focusing wave generation, the local Ursell number at each testing location is calculated by:

$$\mathbf{U}\_{\mathbf{r}} = \frac{\eta\_{\text{max}} \lambda^2}{\mathbf{h}^3} \tag{10}$$

The Ursell numbers at different positions are plotted in Figure 12a,b respectively for the QCWA focusing wave and the QCWS focusing wave. From Figure 12a, it is seen that the Ursell number arrives at the maximum as the wave group approaches the focus position, which implies that the wave nonlinearity plays an important role in the QCWA focusing wave generation. This is also observed

through the extra higher-order harmonics' variation at different positions, as shown in Figure 10. However, from Figure 12b, one can observe that the Ursell number of the QCWS focusing wave keeps rising although the wave group has passed the focus position. It is deduced that the QCWS focusing wave generation is probably dominated by wave components in input frequency range. Additionally, the extra higher-order harmonics have smaller amplitudes in wave spectra of Figure 11 which could be considered as evidence of this deduction.

**Figure 12.** Ursell number changing at different space positions for two generations of focusing waves: (**a**) The QCWA focusing wave and (**b**) the QCWS focusing wave.

#### *3.2. Focusing Wave Generation under Di*ff*erent an and* ω*<sup>n</sup> Parameters*

The QCWA and QCWS focusing waves with different frequency ranges are generated in the physical wave tank. The required maximum displacements (i.e., Smax) of hydraulic piston are summarized in Table 1. For wave groups 1-1 and 2-a, their maximum displacements of hydraulic piston both increase with their lower bonds of wave frequencies being extended, whereas, when the upper bonds of the frequency decrease, for wave groups 2-a, 3-b and 4-c, the corresponding maximum displacements also increase. These comparisons indicate that focusing waves with more low-frequency wave components require wider wavemaker stroke range.


**Table 1.** Required maximum displacements of hydraulic piston.

By comparing the maximum displacements of the QCWA and QCWS focusing waves in each group, we find the QCWS focusing wave requires wider movement range of the wavemaker and this difference becomes more obvious when the frequency range gets wider. For the example of the focusing wave group 1-1 in Table 1, the maximum displacement of the QCWS focusing wave is almost twice that of the QCWA focusing wave. Therefore, the QCWA spectrum is more effective for the generation of focusing waves containing a wider frequency range from the aspect of the wavemaker stroke.

The generated focusing wave crest elevations of the last three experimental groups (i.e., groups 2-a, 3-b, 4-c, in Table 1) are plotted in Figure 13a. The wave crest elevations from the two spectra both raise when their frequency upper bounds are decreased. A similar relationship between the frequency bandwidth and the focusing wave crest elevation is also found in experimental results of the CWA focusing waves [19] and in numerical investigations [24]. We infer that the wave crest elevation is significantly affected by the wave frequency bandwidth but slightly affected by the type of wave amplitude spectrum, as shown in Figure 13a.

Figure 13b compares focusing positions of the QCWA and QCWS spectra. We observed that both waves focus closer to the flap, with the frequency bandwidth narrower. At the same time, the above results also demonstrate that the focusing position changes with the corresponding central frequencies. Li and Liu [50] presented an opposite correlation between focusing position and frequency bandwidth for the focusing waves based on CWA and CWS spectra, but the correlation between the focusing position and the central frequency in their work is similar to that illustrated in Figure 13b. The comparison of the two tested spectra in this study, showed that the focusing position of the QCWA focusing wave is a bit farther from the flap. Therefore, if a focusing wave with a larger wave crest elevation and a farther location from wavemaker is required, the QCWA spectrum is better than the QCWS. In addition, the QCWA focusing wave can also include a wider frequency range for a specific wavemaker with the fixed stroke range, according to Table 1. The QCWA focusing waves will be further investigated in the following sections.

**Figure 13.** Comparison of generating results of the QCWA and QCWS focusing waves: (**a**) Focusing wave crest elevation and (**b**) focusing wave position (a: ω<sup>n</sup> = 1.336 rad/s~3.696 rad/s, b: ω<sup>n</sup> = 1.336 rad/s~3.142 rad/s, and c: ω<sup>n</sup> = 1.336 rad/s~2.094 rad/s).

#### *3.3. Focusing Wave Generation under Di*ff*erent Af and xf Parameters*

In order to investigate the effect of Af and xf on the generating results, focusing wave groups with two different wave frequency ranges are tested in a physical wave tank. In the first experiment group, ω<sup>n</sup> = 0.997 rad/s~3.696 rad/s and Af is chosen as 0.1 m and 0.14 m. For the second group of focusing waves, ω<sup>n</sup> is varied from 1.336 rad/s to 3.142 rad/s with three wave crest elevations, i.e., Af = 0.1 m, 0.14 m, and 0.20 m. All these focusing waves are tested at two different focusing locations, i.e., xf = 30 m and xf = 50 m. The focusing time parameter tf and wave component number N are assigned as 36 s and 32.

Figure 14a,b illustrates tested focusing wave crest elevations under different Af and xf parameters. The generated relative wave crest elevations increase nonlinearly with Af being larger at both specified focusing locations, i.e., the increasing ratio is not a constant. The same conclusion can be also drawn

from the comparison of the CWS and CWA focusing waves by Li and Liu [50] and in experimental investigations on the CWA focusing waves by Baldock et al. [19]. The wave crest elevation increase is generally attributed to the generated nonlinear wave components (i.e., extra higher-frequency components) produced in focusing wave generation. By maintaining Af and varying xf parameter from 30 m to 50 m, it is seen that the tested focusing wave crest elevations show a slight decrease in Figure 14a,b. This may be because the wave energy does not fully focus at the assigned location or a small portion of wave energy is dissipated in focusing wave propagation.

Figure 14c,d demonstrates the actual positions of focusing waves with two wave frequency ranges. As for focusing wave cases with ω<sup>n</sup> = 0.997 rad/s~3.696 rad/s, the larger wave elevation Af results in a forward moving of the focusing point. But the focusing wave group with the other frequency range shows an inversed changing trend. Despite the irregular focusing position variation existing, the overall shift of the focusing position is moderate and the waves approximately focus at 0.9 × xf in the physical wave tank.

**Figure 14.** Generated results of focusing waves under different Af and xf parameters: (**a**) Focusing wave crest elevations for ω<sup>n</sup> = 0.997 rad/s~3.696 rad/s, (**b**) focusing wave crest elevations for ω<sup>n</sup> = 1.336 rad/s~3.142 rad/s, (**c**) focusing wave positions for ω<sup>n</sup> = 0.997 rad/s~3.696 rad/s, and (**d**) focusing wave positions for ω<sup>n</sup> = 1.336 rad/s~3.142 rad/s.

#### *3.4. Focusing Wave Generation under Di*ff*erent xf and tf Parameters*

A series of focusing wave cases are tested against various focusing position and time parameters with a fixed amplitude of Af = 0.1 m. The total number of wave components is 32 and the wave component frequency ω<sup>n</sup> is varied from 1.336 rad/s to 3.696 rad/s. The focusing time parameter is, respectively, set as 25 s, 30 s, and 36 s at each focusing position of xf = 30 m, 40 m, and 50 m.

The measured crest elevations of focusing waves are illustrated in Figure 15a. It is found that when xf = 30 m the wave crest elevations are almost the same for three focusing times and they are about 10% higher than the expected wave crest elevations (i.e., Af 0). As the focusing position parameter is increased to 40 m, the generated wave crest elevations with tf being 30 s and 36 s basically stay the same as those tested at xf = 30 m. However, the wave crest elevation with tf being 25 s shows a significant drop. When the focusing position is further increased to 50 m, the wave crest elevations highly depend on the focusing time parameter and they significantly increase from about 0.6 to 1.0 as the focusing time parameter ranges from 25 s to 36 s. Comparing the wave crest elevations under three different xf parameters, it is found that the wave crest elevations with the focusing time tf being 25 s and 30 s both show a decreasing trend as the focusing position is changed from 30 m to 50 m. Nevertheless, when the focusing time is increased to 36 s, the generated focusing wave crest elevations are largely maintained as the focusing position parameter varies. We deduce that only when the focusing time parameter tf exceeds a threshold, i.e., about 36 s for this testing case, the expected wave crest elevation can be easily achieved in the physical wave tank.

**Figure 15.** Focusing wave generating results under different focusing position and time parameters: (**a**) Focusing wave crest elevations and (**b**) focusing wave positions.

Figure 15b compares the actual focusing positions against theoretical focusing positions under different focusing times. We observed that experimental focusing waves appear in the position of about 0.9 × xf from the flap, as the focusing time tf is 30 s or 36 s. When the focusing time decreases to 25 s, the ratio of actual and assigned positions is below one and it becomes smaller as the given xf increases from 30 m to 50 m. Thus, according to above analysis, a sufficient time length tf for focusing wave generation is suggested to achieve smaller deviation of the actual focusing position.

To find the relationship between the focusing position, xf, and the focusing time, tf, from the wave energy aspect, the required minimum focusing time tcr, i.e., the time for the highest-frequency wave component energy propagating to the specified focusing position, is calculated by the linear wave theory. As for xf being specified to be 30 m, 40 m, and 50 m, the calculated tcr are 22.59 s, 30.12 s, and 37.65 s, respectively. According to the comparison between the tf parameters of the above focusing wave cases and their corresponding tcr results, we deduced that the smaller wave crest elevation and the forward shift of focusing positions, as seen in Figure 15, result from the wave energy not being completely focused at the specified positions. Therefore, for focusing wave generation in a physical wave tank, the assigned time parameter tf should be longer than tcr as predicted by linear wave theory. It should be noted that for focusing wave cases with large crest elevation the time parameter, tcr, predicted by the linear wave theory may not be applicable due to strong wave nonlinearity.

#### **4. Conclusions**

In this paper, two new wave amplitude spectra (i.e., the QCWA and QCWS spectra) used for focusing wave generation in a wave tank have been formulated and experimentally tested. With fixed wave energy, the amplitude and steepness of wave components in these two wave spectra are adjustable. Focusing waves have been generated in the physical wave tank by modifying the wave generator

signal through adding two extensions at the beginning and terminating stages. A series of focusing waves based on the two newly-proposed wave amplitude spectra are tested in physical wave tank, demonstrating that by comparing with the QCWS spectrum the focusing waves generated based on the QCWA spectrum can achieve larger wave crest elevation, and are focused further downstream away from the wavemaker. Additionally, to generate a focusing wave with the same wave parameters, the QCWA spectrum requires a smaller wavemaker movement range. The spectral analysis results show that the wave nonlinearity plays an important role in the QCWA focusing wave generation, whereas the redistribution of wave energy in input frequency range significantly affects the focusing wave generated by the QCWS spectrum. It is also found that two spectra have similar effects on wave crest alternation through adjusting wave frequency range and wave crest elevation parameter, whereas actual focusing positions are highly dependent on the type of the wave spectrum.

Experimental investigations on space-time parameters have also been carried out. We found that the space parameter, xf,, and the time parameter, tf,, jointly affect the generated focusing wave crest elevation and focusing position. The focusing time parameter is constrained by a critical value, which could be predicted by linear wave theory for linear or weakly nonlinear focusing waves. On the basis of this research, future work should be carried out on strongly nonlinear focusing wave generation using these two newly formulated spectra and investigating the effects of the space-time parameters on focusing waves with large-crest elevation.

**Author Contributions:** Investigation, G.X.; validation, H.H.; supervision, Q.M.; writing—original draft; writing— Review and Editing, Q.G.

**Funding:** The National Natural Science Foundation of China (grant No. 51739001 and No. 51909125), the Zhejiang Provincial Natural Science Foundation of China (grant No. LY17A020002), and the Ningbo University Science Foundation (grant No. XYL19030).

**Acknowledgments:** The authors want to thank Fenglai Li and Shizhou You, working in towing wave tank of the Harbin Engineering University. They provided help and support for the focusing wave experiments.

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

#### **Nomenclature**


#### **Appendix A**

Focusing wave cases involved in this research are summarized in Table A1, in which the related parameters and wave amplitude spectrum types of each focusing wave case are listed in detail. For parameters in Table A1, their definitions can be found in the Nomenclature section. From Table A1, it is found that the BFI value of each focusing wave case is quite small, which indicates that the nonlinear self-focusing of wave components almost has little influence on focusing wave generation of this research. The focusing wave generation of this research is dominated by the space-time focusing of wave energy.


**Table A1.** Focusing wave cases in a physical wave tank.

#### **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/).
