<sup>2</sup> Closed boundary condition:

The normal flow rate at the shoreline of the given water body should be 0.

#### 2.2.2. The Euler Model for Residual Current Calculation

As the most important environmental dynamic factor in coastal waters, the residual current plays a crucial role in the transport and diffusion of substances in seawater. In studies on ocean dynamics, Eulerian velocity is usually used to calculate residual currents. The Euler residual current in the ocean can be simply defined as the mean Eulerian velocity, which can be calculated with the following equation:

$$\text{sLI}\_{\text{E}} = \frac{1}{nT} \int\_{t\_0}^{t\_0 + nt} u(\mathbf{x}\_{0\prime} t)dt,\text{ } V\_{\text{E}} = \frac{1}{nT} \int\_{t\_0}^{t\_0 + nt} v(\mathbf{x}\_{0\prime} t)dt\tag{5}$$

where *UE* and *VE* represent the mean Eulerian velocities in directions *x* and *y*, respectively; *n* represents the number of cycles used in the calculation; *t*<sup>0</sup> represents the start time of calculation; *T* represents the current cycle; *u(x*0*,t)* and *v(x*0*, t)* represent the component velocities in directions *x* and *y*. The numerically discrete form of Equation (6) is described below:

$$\mathcal{U}\_{\rm E} = \frac{1}{N} \sum\_{i=1}^{N} u\_{i\prime} \; V\_{\rm E} = \frac{1}{N} \sum\_{i=1}^{N} v\_{i} \tag{6}$$

where *N* = *nT*/Δt and Δt represent the time step of numerical simulation.

#### 2.2.3. Mathematical Model of Water Exchange

The tracer method was used to simulate the degree of water exchange [22–24], where a dissolved non-degradable and conservative substance was set in the sea area, and its concentration diffusion under the action of hydrodynamic force was investigated. For the transport of the tracer, the convection–diffusion equation based on Eulerian substance transport was used, as shown below:

$$\frac{\partial h\mathcal{C}}{\partial t} + \frac{\partial hu\mathcal{C}}{\partial x} + \frac{\partial hv\mathcal{C}}{\partial y} = \frac{\partial}{\partial x}\left(hD\_x\frac{\partial\mathcal{C}}{\partial x}\right) + \frac{\partial}{\partial y}\left(hD\_y\frac{\partial\mathcal{C}}{\partial y}\right) - Fh\mathcal{C} + \mathcal{S} \tag{7}$$

where *C* represents the substance concentration; *Dx* and *Dy* represent the substance diffusion coefficients in directions *x* and *y*, respectively; *F* represents the substance attenuation coefficient, which is zero (*F* = 0) for the conservative substance; *S* represents the point source concentration.

The substance diffusion coefficient was calculated with the following equation:

$$D\_{\mathbf{x}} = \frac{E\_{\mathbf{x}}}{\sigma\_T}; \ D\_{\mathbf{y}} = \frac{E\_{\mathbf{y}}}{\sigma\_T} \tag{8}$$

where *Ex = Ey* represents the horizontal turbulent viscosity coefficient; *σ<sup>T</sup>* represents the Prandtl number, which was determined as 1.0 in this study.

After a certain period of time, the percentage of the total amount of substance diffused from the system to open water divided by the total amount of initial substances in the system should be the water exchange rate of the overall system. The statistical calculation expression is provided below.

$$EX(t\_j) = \left(1 - \frac{\sum\_{i=1}^{N} \mathbb{C}\_i(t\_j) H\_i(t\_j)}{\sum\_{i=1}^{N} \mathbb{C}\_i(t\_0) H\_i(t\_0)}\right) \times 100\% \tag{9}$$

where *EX* represents the water exchange rate; *C* represents the substance concentration; *H* represents the total water depth; *i* represents the node number in the statistical domain; *n* represents the total number of nodes in the statistical domain; *j* represents the time number.

#### *2.3. Grid Creation and Parameter Setting*

#### 2.3.1. Grid Creation

The calculation grid was generated with the Surface Water Model System(SMS 10.1). This grid generation program can realize a flexible and variable resolution in the horizontal direction of the grid and a large gradient, and it can create a highly smooth grid at a location where a flow tends to be generated around an island. In addition, it can partially increase the density in areas with complex terrain, such as coastal areas, estuaries, and wetlands.

The entire two-dimensional simulated domain of Changhai County consists of 32,560 nodes and 63,579 triangular elements. Figure 3a shows the calculation domain and grid distribution of the established two-dimensional model of the sea area near Changhai County. The entire simulated domain of the three-dimensional model [25] consists of 12,539 nodes and 24,281 triangular elements. Figure 3b shows the calculation domain and grid distribution of the three-dimensional model of a small area of Changhai County.

#### 2.3.2. Model Calculation Settings

#### (1) The calculation of hydrodynamic force

The calculation time step of the model was adjusted according to the CFL conditions to ensure that when the model calculation was converged, the minimum time step was 5.0 s. The seabed friction was controlled by the Manning number, with a specific value of 32–42 m1/3/s. In many applications, a constant eddy viscosity can be used for the horizontal stress terms. Alternatively, Smagorinsky proposed to express sub-grid-scale transport by effective eddy viscosity related to a characteristic length scale. The sub-gridscale eddy viscosity is given in [26], and the specific expression is shown below:

$$A = c\_\*^2 l^2 \sqrt{2S\_{ij}S\_{ij}}\tag{10}$$

where *cs* is a constant; *l* is a characteristic length; and the deformation rate is given by

$$S\_{ij} = \frac{1}{2} \left( \frac{\partial u\_i}{\partial x\_j} + \frac{\partial u\_j}{\partial x\_i} \right) \tag{11}$$

The minimum calculation time step of the three-dimensional small-scale local model was 1.0 s. In the vertical grid, the sigma hierarchical function was adopted, and the impact of the rafts on the hydrodynamic force was described with a double-resistance model with the introduction of a secondary drag coefficient, where the frictional resistance of the seabed was controlled by secondary drag coefficient Cf, and the specific expression was determined by assuming a logarithmic profile between the seabed and a point at a distance of *DZb* above the seabed as follows:

$$\mathbf{C}\_{\mathbf{f}} = [\frac{1}{\kappa} \ln \left( \frac{D Z\_{\mathbf{b}}}{Z\_0} \right)]^{-2} \tag{12}$$

where *κ* = 0.4 is the von Kármán constant; *Z*<sup>0</sup> represents the length scale of the roughness of the riverbed; when the boundary surface is rough, *Z*<sup>0</sup> depends on the roughness height, where *Z*<sup>0</sup> *= mks*, the approximate value of *m* is 0.033, *k*<sup>s</sup> is the roughness height, ranging between 0.01 m and 0.30 m, and the value was determined as 0.05 m. In summary, the average value of the secondary drag coefficient was 0.01.

(2) The calculation for the aquaculture area

The location of the selected aquaculture area, the range of the calculation domain, and the distribution of the seabed topography are shown in Figure 1. The aquaculture area is located in a sea area near the Changshan Archipelago in the southeast of Changhai County, and its boundaries are shown in Table 3 below.

**Table 3.** Scope of the aquaculture area.


In the post-aquaculture model, unstructured grids were also used to divide the horizontal calculation domain and locally densify the sea area where the aquaculture area was located, with a grid scale of 30 m. In other areas along the shoreline, the grid scale ranged between 50 m and 100 m; in sea areas far away from the aquaculture area, the maximum grid scale was 400 m, and the calculation domain contained 17,642 triangular grids and 9011 nodes. The Manning field considering aquaculture areas of different densities is shown in Figure 4 below.

**Figure 4.** Distribution of Manning field after the implementation of aquaculture activities.

(3) Assessment of water exchange capacity

In this study, the water exchange rate was used as an index to describe the water exchange capacity of the aquaculture area. A dissolved conservative substance was placed in the sea area where the aquaculture area was located, which would be carried by the water body and could not be degraded. The convection and diffusion of the conservative matter directly reflect the form of movement of the water body. Based on the above considerations, in this study, a conservative substance with a concentration of 1.0 was placed in the aquaculture area; the concentration of substances in open water was set at 0.0; the attenuation coefficient was set as *F* = 0, and the point source concentration was set as *S* = 0. The substance diffusion coefficient was equal to the turbulent viscosity coefficient of water flow (*σ<sup>T</sup>* = 1.0).

#### **3. Results**

#### *3.1. Model Verification Results*

The actual calculation and simulation period of the model was from 0:00, 1 August 2021 to 23:00, 31 October 2021. Figure 5 shows a time-curve-based comparison between the simulated and measured values of the tidal level during the period from 0:00, 8 August 2021 to 23:00, 15 September 2021. Figure 6 shows the comparison between P1 and P2 in terms of flow rate, flow direction, and measured value during the period from 11:20, 13 September 2021 to 12:30, 14 September 2021. To fully represent the calculation results of the numerical model, Figure 7 shows the flow rates and flow direction fields of the top, middle, and bottom layers at the same moment during a spring tide and a neap tide in the three-dimensional calculation model for Changhai County.

**Figure 6.** Verification of flow rates and directions at P1 and P2.

**Figure 7.** Flow field distribution maps of the top, middle, and bottom layers during a spring tide in the three-dimensional calculation domain of a small area of Changhai County.

A comparison between Figures 5 and 6 shows that the results calculated by the model are in good agreement with the measured values, where the error is within an acceptable range. According to the verification results of tidal currents and the flow field distribution maps at different moments, as shown in Figure 5, the mathematical model can reflect the flow field in the sea area near Changhai County in a more realistic manner, indicating that the model has reasonable boundaries and parameters and can be used for the calculation of subsequent working conditions. Figure 6 indicates that the trends of variation in the simulated flow rate and flow direction are generally the same as those of the measured values, where the measured maximum flow rates of P1 and P2 during 13 September and 14 September are 0.606 m/s and 0.491 m/s, respectively, and the simulated maximum flow rates are 0.577 m/s and 0.493 m/s, respectively. In general, the simulated calculated values are in good agreement with the measured values. Figure 7 shows that although grids of different scales were used in the calculation of the two-dimensional and three-dimensional models, the flow rate and flow direction in the entire spatial calculation domain reasonably changed, featuring a strong gradient and no sudden change, which indicates that the model is stable and can be used as the basis for the calculation of subsequent working conditions.

In order to verify the accuracy, the root-mean-square error (R) was used to quantitatively analyze the error between the calculation results of the model and the measured values, where *R* represents the mean deviation between the results of the model and the observed data. *R* is calculated as follows:

$$R = \sqrt{\frac{\sum\_{i=1}^{n} \left(M\_i - O\_i\right)^2}{n}} \tag{13}$$

where *Mi* represents the calculated value of the model; *Oi* represents the observed value; and *n* represents the number of observed values. After calculation, the root-mean-square error in the tidal level of T1 is less than 0.18 m, generally indicating that the simulated value is in good agreement with the measured value, and the model has great accuracy and reliability.

Tracer model verification is indeed a relatively important part of the assessment of water exchange capacity, which has been verified in other studies [27,28], as described below. In these studies, one of the major tasks was the numerical simulation of the transportation of water pollutants in Liaodong Bay, the northernmost bay in the Bohai Sea in China. On

the basis of a comprehensive understanding of the natural conditions of the sea area of Liaodong Bay, the impact of point source input was introduced to establish a convectiondiffusion model for pollutant transport in the sea area of Liaodong Bay, with which the distribution of different nutrient elements in the sea area was simulated, and where major water quality indicators included NH3-N and COD. The accuracy and stability of the model were verified through a comparison between the simulated values and the concentration levels of the elements obtained by field sampling and analysis in the sea area, indicating that the model features a great ability to reproduce and predict the concentration field in the sea area of Liaodong Bay. In respect of the following few years, the distribution of PO4-P, a water quality indicator, in Liaodong Bay, was reproduced, and model verification was performed, providing the simulation results of the hydrodynamic field, half exchange time, and concentration field in Liaodong Bay at different typical moments. Finally, this model was adopted in the simulation and assessment for the identification of marine pollution accidents, and it delivered satisfactory results.

#### *3.2. Simulation Results of Tidal Current and Residual Current in the Aquaculture Area*

The results regarding the distribution of the flow field in the sea area near the aquaculture area are shown in Figures 8–11, indicating that the flow field in the aquaculture area exhibits the characteristics of reciprocating motion, where the main flow direction is NW–SE, and the flow rate magnitude during a flood tide and an ebb tide is 1.0 m/s. In accordance with the simulation result of the hydrodynamic field in the sea area of Changhai County, the Euler residual current field in the sea area near the aquaculture area was obtained. Figure 12 shows the Euler residual current fields in the cycle of a spring tide before and after the implementation of aquaculture activities, indicating that the mean residual current intensity in the aquaculture area was approximately 0.018 m/s; the direction was NE; and there was generally no significant variation before and after the implementation of aquaculture activities.

**Figure 8.** Flow field distribution at high tide.

**Figure 9.** Flow field distribution at low tide.

**Figure 10.** Flow field distribution at the peak of an ebb tide.

**Figure 11.** Flow field distribution at the peak of a flood tide.

**Figure 12.** Residual current distribution in the sea area near the aquaculture area ((**left**): before aquaculture; (**right**): after aquaculture).

In order to quantitatively and clearly reflect the impact of surface roughness (Manning) on the calculation results, the sensitivity of the Manning number was analyzed. For the sea area near the floating raft aquaculture area, which is 180–3650 m from the shoreline, the simulated values of flow rate and flow direction at the peak of a flood tide and at the peak of an ebb tide during a spring tide in one tidal cycle were compared before and after the implementation of aquaculture, where the Manning settings under the two working conditions were as shown in Part 2. The results show that the changes in flow rate and flow direction were generally significant in the study area, especially during the flood tide, wherein the flow rate changed by more than 80% within 750 m in the aquaculture area; the mean change in flow rate was approximately 10%, and the number of points where the flow direction changed by more than 45◦ accounted for around 20% of the total number. This

indicates that, after the establishment of aquaculture activities, as the Manning number of the aquaculture area decreased and the roughness increased, which, together with the effect of flow resistance arising from the raft net, indeed significantly affected the flow rate and flow direction of the sea area near the aquaculture area.

#### *3.3. Simulation Results of Water Exchange Rate*

Based on the aforesaid hydrodynamic simulation results, a mathematical model of water exchange, described in Equations (7)–(9), was used to study the water exchange capacity of the aquaculture area. In addition, in order to enable the model to reflect the variation in water exchange capacity before and after the implementation of aquaculture activities in a clearer and more sensitive manner, the scope of the aquaculture area was appropriately magnified according to the actual sea area, and the calculation domain and grid were rearranged and refined, where the average grid scale was 50 m, and the minimum grid scale for the key areas of interest in the aquaculture area was 30 m.

Figures 13 and 14 show the initial field distribution of the tracer concentration distribution, and the overall water exchange rate–time curves before and after the implementation of aquaculture activities, respectively. Table 4 shows the statistics for the overall water exchange rate–time curves before and after the implementation of aquaculture activities, indicating that the water exchange rate after the implementation of aquaculture decreased compared with that before implementation. Before the implementation of aquaculture, the water exchange rates after 1, 4, and 8 days of water exchange were 27.90%, 61.83%, and 76.48%, respectively; after the implementation of aquaculture, the water exchange rates after 1, 4, and 8 days of water exchange were 22.90%, 53.43%, and 75.23%, respectively.

**Figure 13.** Tracer concentration distribution at the initial time point (the red area is the statistical area for water exchange rate).

**Figure 14.** Water exchange rate–time curves before and after the implementation of aquaculture.


**Table 4.** Statistics for the water exchange rate–time curves before and after the implementation of aquaculture.

#### **4. Discussion**

#### *4.1. Tidal Current Conditions*

The hydrodynamic force calculation results indicate that the hydrodynamic force in the waters near Changhai County is mainly in the NW–SE direction during a spring tide, during which the tidal currents of flood and ebb tides rotate counterclockwise. During the spring tide, the tidal field shows that the tidal direction in the open waters of Changhai County was NW, the tidal field was stable, and the flow rate generally ranged between 0.50 m/s and 0.85 m/s. The nearshore current is a coastal current in essentially the same direction as the shoreline and has a lower flow rate, ranging between 0.2 m/s and 0.4 m/s. This is mainly because the flow rate is significantly reduced by the bottom friction due to the shallow water in the near-shore area. The local coastal waters are affected by the coastline, with a maximum flow rate of 1.2 m/s during a flood tide; during an ebb tide, the flow direction in the open waters is SE, and the flow rate ranges between 0.45 m/s and 0.90 m/s. Affected by the topography, coastal waters generally feature lower flow rates, with a maximum flow rate of approximately 1.1 m/s.

#### *4.2. Conditions of the Aquaculture Area*

The calculation of the characteristics of the flow rate in the aquaculture area shows that flow rates in the aquaculture area usually range between 0.2 m/s and 0.4 m/s. There is a small island in the SE direction in the aquaculture area, so the maximum flow rate in the aquaculture area is up to 0.7 m/s; the flow rate gradually increases from the shoreline to the sea, and it reaches 1.0 cm/s in the part of the aquaculture area that is closer to the shoreline. A comparison with the flow rate before the implementation of aquaculture activities shows that, with regard to the degree of flow resistance imposed by the floating rafts, the variation in flow rate ranges between 2.87% and 84.58% at the peak of a flood tide, and between 2.65% and 20.89% at the peak of an ebb tide from a low-density zone to a high-density zone of the aquaculture area. This indicates that the variation in flow rate caused by the floating rafts in the sea area near the aquaculture area of Changhai County is significantly greater during a flood tide than that during an ebb tide, and the flow resistance rate at the peak of a flood tide is greater than 80%. Therefore, aquaculture operators and marine environmental protection workers should pay attention to the impact of floating rafts for aquaculture. Even in open sea areas, during the setting of the orientation and density of a floating raft aquaculture area, it is crucial to first investigate the hydrodynamic conditions and the impact of aquaculture activities on the hydrodynamic conditions in the sea area, in order to scientifically implement aquaculture activities and rationally determine the layout while protecting the marine environment.

#### *4.3. Residual Current Conditions*

Residual current distribution plays a decisive role in the transport and diffusion of bait, nutritive salts, and other related substances in an aquaculture area. According to a comparison with the residual current before the implementation of aquaculture activities, the extent of variation ranged between 3.01% and 84.74% during a spring tide, and it ranged between 9.46% and 78.50% during a neap tide, indicating that the extents of variation in residual current during tides are essentially the same; they should not be underestimated. Therefore, to accurately understand the distribution of algae and bait in the floating raft aquaculture

area, we must calculate and analyze the residual current in the sea area based on accurate hydrodynamic analysis. In this way, we can understand the characteristics of the transport and diffusion of floating and suspended substances in the sea area in real time, thereby providing guidance for the formulation of aquaculture plans and production activities.

#### *4.4. Water Exchange Conditions*

The quantitative calculation shows that, due to the aquaculture activities, the water exchange rates of the open sea area decreased by 17.92%, 13.59%, and 1.63% compared with those before implementation; moreover, the half-exchange cycle of the water body appeared in 2.3 d and 3.9 d, respectively, before and after the implementation of aquaculture. This indicates that even floating rafts for aquaculture located in an open sea area have a certain impact on the water exchange capacity, and the specific extent of such impact is closely related to various factors, such as the density, size, scope, and location of rafts in the aquaculture area.

#### **5. Conclusions**

In this study, a numerical simulation method was applied to a floating raft aquaculture sea area to quantitatively calculate and assess the changes in the hydrodynamic environment of the open sea area. The model is based on the solution of the three-dimensional incompressible Reynolds-averaged Navier–Stokes equations. Then, the integration of the horizontal momentum equations and the continuity equation over depth for the twodimensional shallow water equations was carried out. In the hydrodynamic model, in order to generalize the impact of rafts on the hydrodynamic force in the aquaculture area, the Manning number of the seabed—namely the seabed roughness—in the two-dimensional mode was changed; in the three-dimensional mode, a double-resistance model of the top and bottom layers was used, with the introduction of a secondary drag coefficient. The final verification and results show that the numerical model proposed in this paper can satisfactorily simulate and predict the hydrodynamic conditions of the sea area near the aquaculture area in Changhai County; the three-dimensional flow field can reflect the variation in the spatially stratified hydrodynamic indexes of the dynamic environment in a more realistic way, which can also reflect the hindering effect of rafts on hydrodynamic force in a more accurate way, and the model features great accuracy and stability. According to the working conditions before and after the implementation of aquaculture activities, the impact of the floating rafts on the hydrodynamic environment and water exchange capacity was compared and analyzed. The results indicate that the flow resistance rate was greater than 80%; the maximum decrease in the water exchange rate was close to 20%. The quantitative results sufficiently show that, even if floating rafts are arranged with a certain density in a completely open sea area, they have a great impact on the hydrodynamic conditions of the sea area. Therefore, aquaculture operators and marine environmental protection workers must pay sufficient attention to the impact of floating rafts for aquaculture on the hydrodynamic conditions of sea area. The establishment of the method in this paper provides a basic model for the rational arrangement of a fully open raft aquaculture area and the scientific determination of breeding density, and it offers a quantitative numerical calculation method for the assessment of the water exchange capacity in aquaculture areas containing flexible objects [29] (such as rafts, vegetation, etc.).It also provides aquaculture operators with technical support in making scientific and effective decisions regarding aquaculture.

In future studies, spatial modeling for floating rafts, mainly including floaters, external aquaculture nets for hanging cages and organisms, will be added, and a fluid–structure interaction-based multiphase flow (volume of fluid, VOF) model will be used to simulate the impact of floating rafts for aquaculture on the dynamic environment and water exchange, in order to provide more accurate and comprehensive technical support for the rational arrangement of aquaculture orientation and the scientific setting of aquaculture density. Furthermore, after the accurate determination of the impact of a raft aquaculture area on

the hydrodynamic conditions, it is also possible for aquaculture operators to reasonably select a site for the installation of bait casting and distribution devices, which can thus help to considerably increase the production efficiency of raft aquaculture, guarantee a stable income for aquaculture operators, and improve the social and economic benefits of raft aquaculture in sea areas in Changhai County and even in other open sea areas with floating raft aquaculture.

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

**Funding:** This research was funded by the Modern Agro-Industry Technology Research System (grant number CARS-49), National Key R&D Program of China (grant number 2018YFD0901604), 2021 Special Program for Marine Economic Development of Liaoning Province, Science and Technology Innovation Fund Program of Dalian (grant number 2021JJ13SN74) and Outstanding Young Sci-Tech Talent Program of Dalian (grant number 2019RJ09).

**Data Availability Statement:** This research did not report any data that are linked to publicly archived datasets analyzed or generated during the study.

**Acknowledgments:** Thanks to the hard work of the field survey researchers, we obtained detailed model verification data. Thanks to the great support and cooperation of the writing and checking researchers, the numerical simulation results in the article could be displayed. The authors would also like to thank all the editors and anonymous reviewers for their helpful comments that greatly improved the quality of the manuscript.

**Conflicts of Interest:** No potential conflicts of interest were reported by the authors. 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**

