The study area is located in the Liaodong Bay, China ( Figure 1A ), one of the main areas for stock enhancement in the Bohai Sea, with a sea area of 18,784.53 km². The Liaodong Bay is situated in the northern part of the Bohai Sea, with an average water depth of 18 m. The bay is influenced by irregular semidiurnal tides, with an average tidal range of 2.7 m. It is one of the most important spawning and breeding grounds for many commercially valuable marine species, such as Portunus trituberculatus and Penaeus orientalis (Xu et al., 2010; Wang et al., 2013). The average annual water temperature in the bay ranges from 8 to 10°C. The sea-ice season typically begins in December and ends in March of the following year (Hou et al., 2020; Wang et al., 2021). The main species for aquaculture release in Liaodong Bay are Portunus trituberculatus, Chinese shrimp and Paralichthys olivaceus. Taking the PTL as an example, the window for stock enhancement and release is from June to July each year. This is because the larvae are prone to freeze and die before this period, while later releases may not be harvested before the onset of the ice season. In addition, the stock enhancement and release is usually carried out during the neap tide period due to the relatively stable marine environment compared to the spring tide period. With this in mind, this study focuses primarily on four neap tide periods occurring in June and July of 2019, namely, early June, late June, early July, and late July.

Portunus trituberculatus is a significant economic crab species, classified under the phylum Crustacea, order Decapoda, and family Portunidae. It is primarily found in the coastal regions of China, Korea, and Japan. In the larval stage, the larvae generally acquire the ability to settle after six moults, and in this process, it takes about 25-30 days for the individuals to become larger and gradually acquire the ability to swim (Chen et al., 2023). Larvae can be released 5 days after hatching, when the cephalothoracic armor width needs to be greater than 6 mm (China, 2014).

Given the importance of the habitat characteristics of PTL in the context of its stock enhancement and release, a numerical model of the marine water environment in the Liaodong Bay was developed in this study. The general model consisted of a three-dimensional hydrodynamic model to reproduce the tidal characteristics of the Liaodong Bay, a three-dimensional temperature and salinity model to simulate the distribution of water temperature and salinity in the Liaodong Bay, and a wave model to reflect the wave propagation characteristics in the Liaodong Bay. The hydrodynamic model was based on the Reynolds-averaged Navier-Stokes equations (RANS) simplified under the hydrostatic assumption. The governing equations of the hydrodynamic model are as follows:

where t is the time; x, y are Cartesian coordinate system coordinates; η is the surface elevation; d is the still water depth; h = η + d is the total water depth; u, v and w are the velocity components in x, y and z directions respectively; f is the Coriolis force coefficient, f = 2ωsinφ (ω is the angular rate of revolution, φ is the local latitude); g is the acceleration of gravity; ρ is the density of water; S_xx , S_xy and S_yx are components of the radiation stress tensor; v t is the vertical turbulent (or eddy) viscosity; p a is the atmospheric pressure; ρ 0 is the reference density of water. S is the magnitude of the discharge due to point sources and ( u s , v s ) is the velocity by which the water is discharged into the ambient water.

The model presented in this study utilized a non-structured triangular mesh ( Figure 1B ), consisting of 14,002 grid elements with a minimum resolution of 100 meters for nearshore areas. The depth and shoreline data used in the model were obtained from the ETOPO1 database provided by the National Oceanic and Atmospheric Administration (NOAA). In addition, the model accounted for four major rivers (the Liao River, the Daliao River, the Daling River, and the Xiaoling River) flowing into the Liaodong Bay, with monthly discharge data provided for each. The predicted tidal level obtained through harmonic analysis was used to provide the open boundary tidal condition.

Apart from the influence of current, water temperature and salinity are also significant factors that cannot be ignored in the growth of PTL (Xu and Liu, 2011; Dai et al., 2014). Based on the validated hydrodynamic model, another model was established to simulate the temperature and salinity distribution characteristics of the Liaodong Bay. The temperature and salinity boundary conditions were provided by monthly average data from the Copernicus Marine Environment Monitoring Center website (https://marine.copernicus.eu/). The model control equations are as follows:

where D v is the vertical turbulent (eddy) diffusion coefficient. H ^ is a source term due to heat exchange with the atmosphere. T s and s s are the temperature and the salinity of the source.

On the basis of temperature and salinity, the food source depletion and flow field heterogeneity caused by wave breaking are potential factors affecting the migration and distribution of marine organisms (Hou et al., 2022). Therefore, the MIKE21-SW model was used to calculate the wave characteristics in the Liaodong Bay, which takes into account meteorological conditions such as wind speed and provides the distribution of the significant wave height in the Liaodong Bay, which is widely used for wave characteristics simulation. The background wind field used meteorological data from the Copernicus Marine Environment Monitoring Officer website (https://marine.copernicus.eu/). The model control equations are as follows:

where N ( x → , σ, θ, t) is the action density, t is the time, x → = (x, y) is the Cartesian co-ordinates; v → = (c_x, c_y, c_θ, c_σ) is the propagation velocity of a wave group in the four-dimensional phase space x → , σ, θ and S is the source term for the energy balance equation.

The energy source term, S, represents the superposition of source functions describing various physical phenomena.

where S in  represents the transmission of energy by wind, S n 1 is the wave energy transfer due to non-linear wave-wave interaction, S ds  is the dissipation of wave energy due to whitecapping, S bot  is the dissipation due to the bottom friction and S surf  is the dissipation of wave energy due to the depth-induced breaking.

Previous studies have shown that fluid shear stress caused by waves and tides affects the transport, attachment and recruitment of marine organism larvae (Shanks et al., 2003; Bolle et al., 2009; Reidenbach et al., 2009; Whitman and Reidenbach, 2012; Hou et al., 2022), while water temperature and salinity directly affect larval feeding and development (Zimmerman and Pechenik, 1991; Nurdiani and Zeng, 2007; McGeady et al., 2021). Existing research and the Chinese technical specifications for the stock enhancement indicate the appropriate growth conditions for PTL (China, 2014; Ge, 2019; Wang et al., 2022), as shown in the Table 1 below.

In the suitability model for the habitat of PTL, it is initially believed that flow velocity, water temperature, salinity, and significant wave height are the main limiting factors for the habitat of PTL. The single-factor Habitat Suitability Index (HSI) was used to quantify the suitability of the target organism in a particular factor. The HSI ranges from 0 (least suitable) to 1 (most suitable), indicating that the factor has no restriction on the survival of the target organism when the HSI index is 1. The Comprehensive Suitability Index (CSI) was calculated based on the HSI of each habitat factor and the main limiting factors for the survival of the organism. The CSI ranges from 0 to 1. In this assessment, it was considered that larvae could settle and survive when the CSI was greater than 0.95. The single-factor habitat suitability index (HSI) and the comprehensive suitability index (CSI) of the larvae on the day of stock enhancement can be expressed as:

where i is the i-th grid cell, t i is the total time in the i-th grid cell that a single suitable condition is met within a day; n is the number of habitat factors that affect the growth of individual organisms; V i , T e m i , S a t i and S W H i are the flow velocity, water temperature, salinity and significant wave height, respectively, which are environmental factors characterizing the unit i. These features are calculated by the habitat environment model. H S I v i , H S I T e m i , H S I S a t i and H S I S W H i are the habitat suitability index based on flow velocity, water temperature, salinity and significant wave height at unit i, respectively.

The Weighted Usable Area (WUA) is calculated by simulating the suitability of larvae to hydrological factors. Specifically, WUA is defined as the summation of the product of the area of each control unit ( a i ) and the CSI within the study region.

During the early stage of stock enhancement and release, the movement of crab larvae is largely determined by ocean currents, due to their small size and weak self-swimming ability (Jones et al., 2015; Giltz et al., 2020). To predict the migration and distribution of the larvae, we developed a particle tracking numerical model. This model considers the larvae as particles (600 in total) that move freely attached to water quality points. After being released within the initial suitable area identified by the habitat suitability model, the transport and diffusion of the larvae follows the following principles:

where a is the drift term; b is the diffusion term; ξ is the random number. To compute the trajectory of Y for a given time discretization, we set the initial value as Y 0 = X 0 and recursively solve for Y using the following equation.

where n=1, 2, 3,...,the value of n is contingent upon the Euler drift coefficient a and the diffusion coefficient b . Δ W = W t − W s ∈N (μ = 0, σ² = Δ _n ) is within the continuous time period ( τ n ≤ t ≤ τ n + 1 ), based on the Gaussian increment of the Wiener process.

We utilized a particle tracking model to assess the migration characteristics of larvae and the likelihood of their persistence within suitable habitats during a specified duration following their release.

The root mean square error (RMSE) is a commonly used evaluation metric for testing the reliability of models.

where S i is the modeling series, O i is the observation series, and N is the total number of data in the series. Given that error measures may not always be indicative of optimal performance for numerical models that simulate natural phenomena such as atmospheric dynamics, ocean circulation, or wave generation and propagation, the concept of the skill model has been introduced to enhance evaluation (Hou et al., 2021).

where M is the modeling series, D is the observation series, D ¯ is the average of the observation data, and N is the total number of data in the series. Performance levels are categorized as: larger than 0.65 excellent; between 0.65 and 0.5 very good; between 0.5 and 0.2 good; smaller than 0.2 poor.

A series of field data were initially collected for the validation of the model reliability, which included tide gauge data from two stations (S1, and S2) obtained from continuous 25-hour measurements at 00:00 (UTC) on May 25, 2022, provided by the China National Marine Information Center; current data from three stations (C1, C2, and C3) measured using a current meter at 02:00 (UTC) on September 18, 2021; monthly averages of seawater temperature and salinity (T1 and T2) for 2021 were obtained from the Copernicus Marine Service website (https://marine.copernicus.eu/); and wave data from one station (W1), obtained from a fixed acoustic wave gauge over a 30-hour continuous observation period at 04:00 (UTC) on May 14, 2016. The measured data collected covered at least one full tidal cycle (semidiurnal tides) to minimize temporal heterogeneity and non-robustness of model predictions.

Based on the comparison between the on-site observation data and the results obtained from model calculations ( Figures 2 , 3 ), it was concluded that the modeled tidal flow values were general agreement with the characteristics of tidal propagation under natural conditions. The simulated results of sea water temperature, salinity, and waves were also found to be in good agreement with the measured data.

The results of the error analysis ( Table 2 ) indicated that the calculated root mean square error (RMSE) of the tidal level has an average value of 0.08 m, while the RMSE of the flow velocity and flow direction has an average value of 0.07 m/s and 15.61°, respectively. The RMSE of the water temperature and salinity are 1.05°C and 1.21‰, respectively. The RMSE of the significant wave height and period are 0.11 m and 0.08 s, respectively. Overall, the observed and simulated values showed good agreement. In addition, the evaluation of skill model for tidal level, flow velocity, temperature, salinity, significant wave height (SWH), and Significant wave period (SWP) were 0.93, 0.90, 0.90, 0.91, 0.88, and 0.90, respectively, indicating excellent performance of the established habitat environment model. Consequently, we concluded that the established numerical model was capable of describing the water environment characteristics (tidal flow, water temperature, salinity, and waves) in the Liaodong Bay.

According to the results of the numerical model of water environment and the habitat suitability model in Liaodong Bay ( Figure 4 ), the habitat suitability index of flow velocity ( H S I v ) in the Liaodong Bay during the four neap tide periods showed a gradual increase from the sea to the land. The flow velocity suitability index of the coastal waters on the north side of Liaodong Bay was generally above 0.95, and other suitable areas were sporadically distributed in shallow waters on both sides of Liaodong Bay, where the flow velocity was relatively low due to the dumpening effect of topographical factors. In the southeastern waters of Liaodong Bay, the flow velocity suitability index was less than 0.20. This was mainly due to the high flow velocity caused by the water depth exceeding 30 meters and the complex coastline. The habitat environment model results showed that the maximum flow velocity in this area during the four periods was 1.87 m/s, which occurred in early June.

Under flow velocity constraints, the maximum extent of suitable habitat for PTL in late June was found to be the largest, approximately 1,488.06 km², accounting for 7.92% of the total area of Liaodong Bay. Conversely, the optimal area was found to be the smallest in early June, at approximately 1,337.49 km², accounting for 7.12% of the total area. The suitable area in early July was 1,371.96 km² and in late July was 1,472.04 km², and the flow rate suitable area accounted for about 7.30% and 7.84% of the whole Liaodong Bay area, respectively. It is evident that under the influence of tidal current stress, the habitat in late June is relatively more conducive to the growth of PTL.

The computational model results showed temporal variations in water temperature within the Liaodong Bay. Water temperature exhibited a significant increase from June to July. At the spatial scale, water temperature gradually decreased from the land towards the ocean. The suitability analysis results ( Figure 5 ) indicated that there were unfavorable temperatures for the growth of PTL in the early June and late July, while the habitat suitability index of water temperature ( H S I T e m ) for the entire Liaodong Bay was 1.0 in late June and early July, indicating that water temperature was highly suitable for the growth of PTL in late June and early July.

In early June, the water temperature in the southern waters of the Liaodong Bay was below 16°C, which limited the growth of the PTL, resulting in H S I T e m of less than 0.1 in the region. However, in late July, the water temperature in the nearshore waters of the northern part of the Liaodong Bay exceeded 28°C, unconducive to the growth of the PTL. As a result, the suitability index for this area was lower. Under the constraint of water temperature, in early June, the suitable growth area for the PTL was 11,873.52 km², accounting for 63.21% of the entire Liaodong Bay area. In late July, the suitable growth area for the PTL was 18,591.48 km², accounting for 98.97% of the entire Liaodong Bay area. In late June and early July, the proportion of suitable area reached 100%. Therefore, it can be seen that under high temperatures, the entire Liaodong Bay is suitable for the growth of the PTL in late June and early July.

According to the distribution of the suitability index for salinity ( H S I S a t ), the majority of the areas in the Liaodong Bay were highly suitable for the growth of PTL. The areas with lower suitability were mainly located in the coastal waters on the north and south sides of the Liaodong Bay. In the nearshore areas on the northern side of the Liaodong Bay, the influx of freshwater from rivers resulted in localized areas of water with a salinity of less than 20‰. On a temporal scale, the H S I S a t was inversely proportional to the river discharge ( Table 3 ). In early June, the relatively low river discharge led to relatively high salinity levels in the estuarine area. In early June, the relatively small river discharge resulted in relatively higher salinity in the estuarine area. However, with the onset of the rainy season, the increasing river discharge gradually reduced the salinity in the nearshore estuarine area, with the maximum river discharge in late July resulting in the lowest H S I S a t in the northern nearshore area of Liaodong Bay. In contrast, on the southern side of Liaodong Bay, due to the deeper water depth, the salinity exceeds 32‰, which limited the growth of PTL.

Under salinity influences, the optimal growth area of PTL was largest in late June ( Figure 6 ), with an area of approximately 18775.53 km², accounting for 99.95% of the entire Liaodong Bay. In early June, the optimal growth area was approximately 18,740.61 km², while in early July, it was approximately 18,772.20 km², and in late July, it was approximately 18,646.65 km², accounting for 99.77%, 99.93%, and 99.27% of the entire Liaodong Bay, respectively. Overall, the salinity of water in Liaodong Bay had little effect on the growth of PTL in June and July.

When considering the growth restriction of PTL by waves alone, we found that the suitability index for waves ( H S I S W H ) in the Liaodong Bay was 1 in all four periods, meaning that the significant wave height (SWH) was less than 0.5 m. This indicates that wave conditions in the entire Liaodong Bay are suitable for larval growth during these four periods. Therefore, we regarded waves as a non-primary limiting factor when calculating the comprehensive suitability index (CSI).

The results of the habitat suitability model ( Figure 7 ) showed that under the combined effects of flow velocity, water temperature, salinity, the suitable areas for the stock enhancement of PTL were mainly distributed in the north of Liaodong Bay. In both June and July, the C S I of PTL showed a trend of gradually increase from sea to land. The C S I of nearshore areas generally exceeded 0.85, indicating that the milder natural conditions in nearshore areas were adequate for the survival of PTL. Based on the results of the C S I , the suitable growth area (CSI>0.95) for the PTL was more concentrated in late June, mainly in the northern waters of the Liaodong Bay.

In terms of suitable PTL growth area, it was largest in late June, reaching 3,217.77 km², while the suitable areas in early June, early July, and late July were 2,503.89 km², 2,253.51 km², and 2,121.73 km², respectively. Therefore, it can be inferred, based on the consideration of the operational space for stock enhancement, the best period for stock enhancement is during the neap tide of late June because it presents the largest suitable area for such a purpose. Subsequently, we initially identified suitable area by determining those with a C S I greater than 0.95 in late June.

After determining the initial suitable area, our interest was in ascertaining the true survival rate of the larvae placed post-release. This was because the distribution of suitable areas in real environments is subject to dynamic changes in environmental factors. Thus, it is crucial to determine whether the larvae will detach from the suitable area during the process of migration along with currents after being released.

Within the initial suitable area ( Figures 7B , 8A identified by the habitat suitability model, larvae (600 in total) were released at even intervals of 2400 meters. The migration trajectory model subsequently predicted the migration positions of these larvae on the 10th, 20th, and 30th days ( Figures 8B-D ). It should be noted that each larva represents the migratory characteristics within a range of 2400m x 2400m around its initial location. The results indicate that 10 days after release, the suitable area decreased, with 26.35% of PTL modelled to be outside the suitable area. The overall depth of larvae decreased compared to the initial time, with an average depth of about -2.12 m (below the surface). After 20 days of release, PTL tended to aggregate in the nearshore area, with an average depth of about -0.89 m; approximately 32.67% of the larvae were located outside the suitable area at this time. When the larvae were released for 30 days, with the reduction of the suitable area, more larvae were detached from it, accounting for approximately 46.81% of the total number of larvae.

The observation implies that despite having selected regions with superior environmental suitability for stock enhancement and initial release, dynamic fluctuations in environmental conditions may cause the PTL to struggle to survive if they move/are moved outside the favorable area. Therefore, our initial estimation of the suitable area may have been overestimated. In the light of this, we evaluated the actual available area for stock enhancement of PTL.

In reality, the survival probability of the PTL is only high when they stay within the suitable habitat for a sufficiently long period of time. Specifically, the probability of larvae escaping from the suitable area after initial release indicates whether the surrounding area of 2400m×2400m around the initial location of the larvae is truly suitable for stock enhancement. If the larvae stay in the suitable area for more than 95% of the time within 30 days after release, the area of 5.76 km² (2400m×2400m) represented by the initial release location is considered as the true available area for natural proliferation of larvae and release of additionally reared ones.

The larvae are released within the initially designated suitable area five days after hatching, with settlement typically requiring 25 to 30 days. We assumed that larvae do not have the ability to swim autonomously prior to settlement in this study. For the purpose of quantification, we tabulated the total initial areas that can be utilized for larval stock enhancement when the proportion of time that larvae spend within the suitable area during a 30-day period reaches 50%, 75%, 95%, and 100%, respectively ( Figure 9 ). The results indicate that a total of 424 larvae satisfied the condition of spending 50% of their time in the suitable area, corresponding to a viable area for growth and release of approximately 2,442.24 km², which represents 75.89% of the initial suitable area and 13.00% of the entire Liaodong Bay area. When the percentage of time spent in the suitable area was 75%, the available area for stock enhancement decreased to about 2004.48 km², accounting for 62.29% of the initial suitable area and 10.67% of the entire Liaodong Bay. Nevertheless, when the percentage of time spent in the suitable area reached 95%, only 33.67% of larvae (202 individuals) met the conditions, and the actual available area for stock enhancement decreased to 1163.52 km², which accounted for 36.16% of the initial suitable area and 6.19% of the entire Liaodong Bay. After being released at their initial area, 147 larvae remained within the suitable area throughout the entire observation period. The actual available area for stock enhancement was only 846.72 km², accounting for 26.31% of the initial suitable area and 4.51% of the entire Liaodong Bay.