Assessment of bivalve carrying capacities and seeding densities in aquaculture areas of Jiaozhou Bay, China, using ecological modeling and the food balance

The increasing growth of the bivalve-culture industry necessitates a scientific assessment of the carrying capacity (CC) of this aquaculture system, in order to ensure sustainability. Here, using a marine plankton ecological model and the food balance, bivalve CCs were estimated for aquaculture zones in Jiaozhou Bay (JZB), an important aquaculture center in China. The results indicated that CCs mainly depended on bivalve filtration rates, which were in turn regulated by temperature. The CCs of Manila clams, razor clams, and scallops were lowest in August – September, while the CCs of oysters and mussels were lowest in October. For all taxa, CCs were higher in areas with higher primary production. Further evaluation indicated that present cultures were saturated, with high mortality rates. To improve culture quality and efficiency, the seeding density of the Manila clam, the most widely farmed species in the JZB, should be reduced to (cid:1) zones in the JZB. Our results indicated that CCs were lower in the summer and autumn because of high FRs and high temperatures; areas with high primary production had higher CCs; and bivalve stocks in the JZB exceeded maximal CC. We recommend reduc-tions in average seeding densities to (cid:1) 700 ind./m 2 for Manila clams, 345 ind./m 2 for razor clams, 60 ind./m 2 for oysters, 165 ind./m 2 for mussels, and 280 ind./cage for scallops. Our study provides references for future estimations of aquaculture CCs and reasonable seeding densities in other bays.


| INTRODUCTION
In recent decades, Chinese mariculture has developed rapidly. The area devoted to mariculture increased from $4.5 Â 10 5 ha (1 ha = 10 4 m 2 ) in 1991 to $2.1 Â 10 6 ha in 2010, while Chinese production increased from 1.9 Â 10 6 t to 1.5 Â 10 7 t over the same period (Jiang, Mu, & Yao, 2013). At present, 84.5% of all mariculture production worldwide originates from China (Food and Agriculture Organization, 2018). However, intensive aquaculture has had a negative impact on the ecological environments of many Chinese bays, leading to increased mortality rates, intensified disease outbreaks, and reduced condition among the cultured organisms (Ma et al., 2013). To maintain aquaculture health and sustainability, aquaculture scale and density must be within the carrying capacity (CC) of a given ecosystem. CC is defined as the yield or density at which production is maximized without negative ecological impacts (Duarte et al., 2003).
Previous studies have calculated bivalve CCs using several different methods. First, some authors have determined CC based on historical statistics. That is, culture density is equivalent to CC when production is maximized; above this CC, bivalve growth decreases, and mortality increases (Filgueira, Comeau, Guyondet, Mckindsey, & Byron, 2015). Second, bivalve CC has been estimated based on hydrographic conditions, food availability, and feeding characteristics (Garen, Robert, & Bougrier, 2004). Third, bivalve CC has been determined based on a material equilibrium, in which food supply is exactly sufficient for bivalve growth (Byron, Link, Costa-Pierce, & Bengtson, 2011). Fourth, bivalve CC has been computed using an ecological modeling approach, integrating the biochemical processes of nutrients, plankton, and bivalves (Guyondet, Roy, Koutitonsky, Grant, & Tita, 2010). The first three of these methods depend on measured data, which may be limited, and they often lack timeliness, which inhibits the assessment of dynamic CCs that vary temporally (Filgueira et al., 2015;Ibarra, Fennel, & Cullen, 2014). Use of the fourth method poses different challenges, including defining acceptable ecological impacts, identifying tipping points associated with ecological resilience, and parameterizing bivalve biological processes (Ibarra et al., 2014). For an aquaculture area, determining a reasonable method is the prerequisite for assessment of CC.
The culturing stocks in Jiaozhou Bay (JZB), an important center of bivalve aquaculture in China, have become saturated as culture scale and seeding density have increased, which caused poor growth and heavy mortality (N. Yu, Wang, Lu, Li, & Zhang, 2013). As a result, bivalve production has decreased or even halted (N. Yu et al., 2013;M. Zhang, 2008). For an example, under a seedling density of 2,500 ind./m 2 , the mortality rate for 2-3-year-old Manila clams reached 47% because of deficiency in food and oxygen (M. Zhang, 2008). In this paper, we estimated bivalve CCs of aquaculture zones (HD, JZ, HW, HS, and HE) in JZB by using the marine plankton ecological model in conjunction with food-balance relationships. This joint approach avoids the data shortages associated with assessments based purely on the material equilibrium, and it also overcomes the difficulties associated with assessments based purely on models. Firstly, the ecological model was used to determine spatiotemporal biomass and primary production of aquaculture zones. Then, these model results were used to compute dynamic CCs of bivalves (Manila clams, razor clams, scallops, Pacific oysters, and mussels) based on the food-balance relationships. At last, seedling densities of these bivalves were investigated with respect to bivalve mortality rates and the estimated CCs, in order to identify optimal densities and ultimately make recommendations.

| The ecological model
We selected model variables according to the principle that key biochemical processes should be reflected, but the model should not be too complicated (Chen, 2003;X. Liu, Li, & Yuan, 2010). First, for the purposes of estimating bivalve CCs based on the food-balance relationship, the ultimate bivalve food source was considered equivalent to primary production. Second, because detritus is not the primary food source for bivalves, detritus was not measured in our survey nor considered in our model. Third, because our aim was to estimate CCs rather than residual capacities by using the modeled biomass and primary production in aquaculture zones, the model did not consider the biological processes of the bivalves. Thus, we adopted a nutrient-phytoplankton-zooplankton (NPZ) model. The biochemical processes and parameterization of the model are shown in Appendix A.

| Model configuration and data usage
Hydrodynamic modeling was performed using the Princeton Ocean Model (Mellor, 2004). Only tidal forces were considered because tidal currents are the prominent physical processes in the JZB (X. Liu et al., 2010). The model domain covered the entire bay ( Figure 1). The open boundary was at the bay mouth, where the water-level forcing was calculated as follows: where A i , θ i , and ϕ i represented the amplitude, phase lag, and initial phase, respectively, of M 2 , S 2 , O 1 , and K 1 component tides. These values were obtained based on harmonic analysis (Huang & Huang, 2007)   Terrestrial nutrients enter the JZB via six rivers ( Figure 1); the associated nutrient input data were obtained from K. Li, Zhang, Li, Zhang, and Wang (2015).

| The food-balance relationship used to estimate CCs
The primary bivalve species cultured in the JZB are Manila clams, Ruditapes philippinarum; razor clams, Sinonovacula constricta; scallops, Chlamys farreri; Pacific oysters, Crassostrea gigas and mussels (Mytilus edulis). Of these, the two clams are bottom cultured, while the oysters, scallops, and mussels are suspension cultured. Our survey divided the aquaculture area into five zones: HE, HS, HW, JZ, and HD ( Figure 1). All five species are cultivated in zones HE, HS, and HW; only Manila clams, razor clams, and oysters are cultivated in zone JZ; and only Manila clams and razor clams are cultivated in zone HD.
F I G U R E 1 Configurations of the model domain, observation sites, benthic sampling sites, and aquaculture zones. HE, HS, HW, JZ, and HD represent Hongdao East, Hongdao South, Hongdao West, Jiaozhou, and Huangdao culture zones, respectively; DX and DY are benthic sampled sites With a known feeding rate and food supply, the maximum number of bivalves supported by the food supply (i.e., the CC) at any given time can be estimated based on the food balance, which assumes an equilibrium between food consumption and supply (Byron et al., 2011). Given primary production PP (representing food supply, unit: mgC/[m 2 hr]) and ingestion rate IR of individual bivalve, the CC, represented with the number of individuals per area, in a given zone was calculated as CC ¼ PP=IR where IR = FRÁkÁP, P was the algal concentration, k (=45) was the mass ratio of carbon to Chl-a in algae, and filtration rate (FR) was the volume of water filtered by the bivalve per unit time (Ehrich & Harris, 2015).
For the suspension-cultured bivalves, PP was equal to algal primary production (PP1). The bottom-cultured bivalves, however, are supplied with food by both algae and microphytobenthos. Therefore, the PP of these bivalves was PP1 plus benthic primary production (PP2).
With the above formula, modeled results (primary production and algal concentration), and observed data (benthic primary production and filtration rates), we computed CCs of aquaculture zones for five species in different sizes in 1 year. Based on the estimated CCs and the mortality rates of bivalves, the stocking densities were evaluated so as to present the suggested seeding densities.

| Observations of benthic primary production
By using 14 C isotope tracer method according the specification for oceanographic survey (SAC, 2008), benthic primary production was measured in seabed samples (0-1 cm columns) taken at sites DX and DY ( Figure 1) in February, April, August, and October in 2000. Benthic primary production in February, April, August, and October was 4.06, 19.77, 8.33, and 47.88 mgC/(m 2 hr), respectively, at site DX, and was 4.45, 9.87, 7.02, and 8.65 mgC/(m 2 hr), respectively, at site DY. Monthly production was estimated by the linear interpolation of the observations, and the contribution of an observed site to a zone was estimated based on the inverse proportion of the distance between the site and the zone.

| FRs of Manila clams
Manila clams are classified into size groups based on size: shell lengths of 2-2.6, 2.6-3.3, and 3.3-4 cm (J. Zhang, Fang, Sun, & Zhao, 2005). The clearance rate (CR), which corresponds to the volume of water filtered per unit time in which food particles are completely ingested, has been studied in Manila clams. J. Zhang et al. (2005) reported the CRs of Manila clams at 21 C, while Dong, Xue, and Li (2000) investigated the effects of temperature on CR. J. Wang et al. (2006) reported the FRs of Manila clams at 15.4 C (seen in Table S1). By extending the results of Dong et al. (2000) (Table S1) in a linear fashion, we estimated that CRs at 4 and 28 C were 0.03 and 0.58 L/hr, respectively. Then, the average CR at 21 C and FR at 15.4 C were obtained. Assuming that the effects of temperature on CR and FR are identical, we combined the results of J. Wang et al. (2006) and Dong et al. (2000) to determine the average FR at 21 C and the ratio of average FR to average CR at 21 C. The ratio was multiplied by CR to yield the FR at 21 C. Finally, FRs at various temperatures were obtained using linear interpolation based on the effects of temperature on FR.
The FRs were linearly interpolated to predict FRs at other temperatures. Fang et al. (1996) divided scallops into three groups based on shell size (3-4, 4-5, and 5-6 cm) and measured the FRs of scallops from Sungo Bay, China, a bay near the JZB. We estimated FRs for the scallops in the JZB taking into account the temperature difference between the JZB and Sungo Bay (Table S3). Kuang et al. (1996) classified Pacific oysters into six size groups based on dry weight (0.12, 0.21, 0.42, 0.66, 1.1, and 1.39 g) and measured the FRs at different temperatures (Table S4). We used linear interpolation based on these results to estimate FRs for different size oysters at different temperatures.

| FRs of oysters
2.9 | FRs of mussels J. Wang et al. (2006) classified mussels into three size groups based on shell length (3.9, 4.8, and 6.5 cm) and measured FRs for three sizes at $17.79 ± 0.87 C (Table S5). However, as FRs were not measured at different temperatures, we were unable to estimate the effects of temperature on mussels based on this previous study. However, as mussels and oysters have similar habits, we assumed that the relationship between oyster FR and temperature was identical to that between mussel FR and temperature. We then estimated mussel FRs at various temperatures by the linear-interpolation of the measured values.
3 | RESULTS Figure 2 shows the modeled and observed monthly Chl-a levels at each site, while Figure 3 shows the monthly variations in average temperature in the bay and the average primary productions in each zone. From winter to spring, as temperature and solar radiation increased, phytoplankton biomass gradually increased. After May, Chl-a levels and primary production decreased slightly as algal metabolism and ingestion by zooplankton increased. During the summer, Chl-a levels and primary production increased continuously because of favorable nutrient, temperature, and solar radiation conditions. Phytoplankton biomass and primary production decreased in autumn. In all seasons, Chl-a levels in the northern area of the bay were higher than those in the central part of the bay and at the bay mouth.

| Simulated phytoplankton biomass
The modeled Chl-a levels were generally consistent with observed values in the bay. However, annual average observed Chl-a (4.08 mg/m 3 ) was higher than that predicted by the model (3.35 mg/m 3 ). This discrepancy was mainly attributed to the maximums recorded in February. Excluding the values measured in February, the annual average observed Chl-a (3.3 mg/m 3 ) was consistent with the modeled value (3.5 mg/m 3 ). In addition, the maximum Chl-a levels in the bay were predicted by the model in August (in the south) and September (in the north). In contrast, excluding February, measured Chl-a levels were highest in June at Sites 3 and 11, in October at Sites 8-10, and in August at all other sites. Overall, measured Chl-a levels were highest in August-October.
Deviations between measured and simulated data are common, and may be because of a number of factors, including survey time and test method. Notably, at northern sites, observed Chl-a levels were highest in February (while quite low in January and March), even significantly greater than those in the summer. Wu, Sun, Zhang, and Zhang (2004) and Sun et al. (2005) reported that, in February, Chl-a peaked, but primary production remained low.
To date, the underlying biological causes of these observations have not been clarified.

| Manila clam CCs
Large monthly variations in CCs of Manila clams were estimated (Figure 4). The CC of large clams was lower than that of small clams, but the CCs of clams of different sizes varied similarly in response to seasonal changes in temperature. At higher temperatures, FRs were higher, and CC was lower. CCs for all sizes of Manila clams were higher in the winter and lower in the summer. In winter, maximum CC for clams of 2-2.6, 2.6-3.3, and 3.3-4 cm were 4,534, 3,364, and 2,744 individuals (ind.)/m 2 , respectively, in the JZB. As temperatures increased in the spring, CC continued to decline, even though primary production increased, because clam metabolism improved significantly. CC was lowest from August to September, with CCs of 630, 468, and 382 ind./m 2 for clams of 2-2.6, 2.6-3.3, and 3.3-4 cm, respectively. As temperatures decreased in autumn, CC gradually increased. It is thus clear that the dominant environmental factor affecting CC was temperature.
Throughout 1 year, the highest is $7 times greater than the minimum. Figure 4 shows the CCs of five culture zones in May-November (CCs in the other months were not presented because they were significantly high and not used in assessing stocking densities). In zones HE, HS, HW, JZ, and HD, the minimum CC for 2-2.6 cm clams was 669, 657, 535, 628, and 744 ind./m 2 , respectively; the minimum CC for 2.6-3.3 cm clams was 497, 488, 396, 466, and 552 ind./m 2 , respectively; and the minimum CC for 3.3-4 cm clams was 405, 398, 323, 380, and 450 ind./m 2 , respectively. In the summer, zone HD had the highest CCs and primary production, while zone HW had the lowest CCs and primary production.

| Razor clam CCs
Because razor clam CCs are highest in the winter, CCs between May and November are shown ( Figure 5). CC varied significantly among months and clams of different sizes. However, in general, larger clams had lower CCs. Seasonal variations in CC followed similar patterns in clams of different sizes. Because filtration is affected by temperature, CCs decreased gradually through the spring and summer, and CCs increased gradually through the autumn and winter. In zone JZ, CCs were lowest in August-September: for 3.9, 5.4, and 6.3 cm

| Oyster CCs
Oyster CCs varied with size and month. Because larger oysters have lower FRs, the CC of larger oysters is lower than the CC of smaller oysters at any given time point. Seasonal variations in CC were consistent among oysters of different sizes: CCs were higher in the winter and lower in the summer and autumn. Figure 7 shows CCs between May and November, when CCs were significantly lower than those in other months. In zone HE, CCs were lowest in

| Mussel CCs
Mussel CCs varied with size and month. In general, bigger mussels had lower FRs and lower CCs. Seasonal variations in CC were consistent among mussels of different sizes: CCs were higher in the winter and lower in the summer and autumn. The maximum CC was observed in February and was $11 times greater than the minimum

Month Month
Carrying capacity (ind/m 2 ) F I G U R E 5 Carrying capacities of razor clams of different sizes CC. Figure 8 shows the CCs of mussels of different sizes between May and November. Similar to oysters, mussel CC was lowest in zone HE in October, but lowest in September in all other zones. Minimum CCs for 3.9, 4.8, and 6.5 cm mussels were 142, 115, and 108 ind./m 2 , respectively, in zone HE; 157, 128, and 120 ind./m 2 , respectively, in zone HS; and 126, 103, and 96 ind./m 2 , respectively, in zone HW. Manila clams are harvested, the maximum production in zones HE, HS, HW, JZ, and HD is predicted to be 29,820, 29,280, 23,760, 27,960, and 33,120 kg/ha, respectively (an average of 28,080 kg/ha). If 3-year-old Manila clams are harvested instead, the maximum production in zones HE, HS, HW, JZ, and HD is predicted to be 40,500, 39,800, 32,300, 38,000, and 45,000 kg/ha, respectively (an average of 38,200 kg/ha). Given the 9,600 ha in the JZB devoted to Manila clam cultivation (Ren et al., 2006), the maximum production of 2-and 3-year-old clams is estimated to be 269,570 and 366,720 t/year, respectively. In recent years, Manila clam production in the JZB has averaged 33,750 kg/ha (a total of 324,970 t) (M. Zhang, 2008), which is in the range for 2-3-year-old clams calculated above.
Therefore, if only 3-year-old clams were harvested, output and income might both increase. That is, the selective harvest of 3-year-old clams might both reduce culture costs and increase economic benefits.
In the JZB, Manila clams are seeded at density of 2,500 ind./m 2 ; the mortality rate for 1-2-year-old clams is 35%, and that of 2-3-year-old clams is 47% (Guo, Ren, & Yang, 2005). Based these data, 411 ind./m 2 remained 2 years after sowing, which was less than the CC of 2-year-old clams (468 ind./m 2 ) but more than the CC of 3-yearold clams (382 ind./m 2 ). It is thus evident that seeding clams at a density 2,500 ind./m 2 does not significantly improve output, but instead increases mortality while decreasing individual quality and fatness, ultimately increasing culture costs. Thus, seeding densities should be appropriately reduced. According to expected survival rates (80% at 2-years-old and 70% at 3-years-old; Guo et al., 2005), we recommend seeding densities (shown in Table 1) in zones HE,HS,HW,JZ,and HD of 723,710,576,678,and 804 ind./m 2 , respectively (an average of $700 ind./m 2 ).
We predicted that razor clam mortality rates were similar to those of Manila clams: 20% at 2-years-old. Thus, we recommend a razor clam seeding density (shown in Table 1)

| Culture densities of scallops
A scallop cage (5.5 m 2 ) has 8 layers with $120 ind./layer; thus, scallop culture density in the JZB is 960 ind./cage (mortality >50%) and the total output is 5,782 t. However, based on the estimated CC of 5-6 cm scallops, we F I G U R E 8 Carrying capacities of mussels of different sizes predicted a yield of 4,883 t in the JZB. Because scallop density exceeds the supported CC, mortality is high. Therefore, scallop culture density should be reduced. M. Zhang (2008) reported that scallop mortality was only 5% if seeded at a density of ≤50 ind./layer, and that mortality increased significantly above this threshold density. C. Li et al. (2011) reported that the seedling density in the JZB and nearby bays was $30 ind./layer in the early 1980s and increased to 50 or even >100 ind./layer in the 1990s, leading to a decrease in algal abundance: algal abundance was 1,000-2,000 Â 10 4 cell/m 3 in 1988, but only $400 Â 10 4 cell/m 3 after 1995. In 2006, we collected an average of 405 Â 10 4 phytoplankton cells/m 3 in the JZB. Overstocking leads to decreased individual quality and high mortality.
In some years, mortality rates reached 90% (C. Li et al., 2011). Y. Wang, Li, Qiu, Zhang, and Wang (2001) reported that mortality rates in these bays were 50-95% between 1995 and 2001, with deaths mainly occurring in the summer. Accordingly, Y.  recommended seeding densities of <40-50 ind./layer, and showed that when seeding density in nearby Sungo Bay decreased to 20 ind./layer, both the ecological environment and the culture efficiency were improved. R. Yu, Lin, and Bi (1994) and Z. Yu et al. (2010) recommended a seeding density of 25-30 ind./layer in a cage.
In the aquaculture zone of the JZB, scallop spats are placed in cages in March and harvested in December as 5-6 cm adults. Therefore, 4-5 cm scallops (September growth) were chosen to determine reasonable seeding densities.
The estimated CC of 4-5 cm scallops in zones HE, HS, and HW was 275, 286, and 226 ind./cage, respectively, equivalent to 34, 36, and 28 ind./layer, respectively. Given a 5% mortality rate, the suggested seedling density (shown in  (Nunes et al., 2003). Gong (1995) reported that the temporary survival rate of seedlings was 75%. Based on these data, we calculated that the survival rate from seeding to harvest was 58%. Thus, the recommended seeding densities (shown in Table 1) in zones HE, HS, HW, and JZ are 62, 68, 55, and 59 ind./m 2 , respectively. The average recommended density across the bay is thus 59 ind./m 2 , consistent with that reported by Nunes et al. (2003) in Sungo Bay.

| Culture densities of mussels
In the JZB, mussels are seeded around November and are cultured more than 1 year; the best harvest time is January-March. Our survey indicated that mussels are usually harvested after reaching a shell length of $6 cm because mussel farmers emphasize economic value and fatness. As CC was lowest in September-October, when mussels are growing rapidly but are not yet ready for harvest, we used the CCs of 4.8 cm mussels to calculate optimal seeding densities. The calculated densities were 115, 128, and 103 ind./m 2 for zones HE, HS, and HW, respectively (115 ind./m 2 on average). Q. Liu (1993) reported that high seeding densities may lead to mortality rates of 60-70% and may significantly reduce individual size and weight; it was suggested that the seeding density should be limited to ≤166 ind./m 2 . Given that 30% of all mussels may die between seeding time and harvest time, we recommend mussel seeding densities (shown in Table 1)

| Culturing arrangement
To maintain the ecological balance and promote healthy aquaculture, rational arrangement in aquaculture systems is prerequisite, so as to increase organismal reproduction and enrich the aquatic food chain (Powers, Peterson, Summerson, & Powers, 2007). In most of the aquaculture zones of the JZB, the water is shallow and the sediment type is of a sandy one, which is advantageous to benthic primary production. This good condition satisfies bottom culture in a large scale with relatively high density. At present, Manila clams occupy 93.8% of the total aquaculture area and represent 95.8% of the total production. This arrangement generally accorded with the condition. Only in the deeper water, scallop, oyster, and mussel are suspended cultured. Besides, culture structures may be adjusted by introducing polyculture (such as kelp and bivalve), choice seafood culture, or fish-cage culture.

| CONCLUSIONS
Here, we presented a method of CC estimation that can be used to assess the sustainability of a mariculture system.
Our approach reflected spatiotemporal variations in bivalve CCs by combining an ecological model and food-balance relationships. We estimated bivalve CCs and evaluated seeding densities in aquaculture zones in the JZB. Our results indicated that CCs were lower in the summer and autumn because of high FRs and high temperatures; areas with high primary production had higher CCs; and bivalve stocks in the JZB exceeded maximal CC. We recommend reductions in average seeding densities to $700 ind./m 2 for Manila clams, 345 ind./m 2 for razor clams, 60 ind./m 2 for oysters, 165 ind./m 2 for mussels, and 280 ind./cage for scallops. Our study provides references for future estimations of aquaculture CCs and reasonable seeding densities in other bays.

SUPPORTING INFORMATION
Additional supporting information may be found online in the Supporting Information section at the end of this article.