Analysing ecological carrying capacity of bivalve aquaculture within the Yellow River Estuary ecoregion through mass-balance modelling

: As the largest aquaculture producer in the world, China is facing the challenge of maintaining sustainability while continuing to develop the aquaculture industry to meet socio-economic needs. Models of trophic structure and energy flow can be used to analyse ecological carrying capacity in order to determine whether a large and rapidly increasing aquaculture industry potentially puts sustainable development at risk. The Yellow River Estuary ecoregion in Shandong Province, China, is an ecologically important region, with extensive bivalve aquaculture that is increasing rapidly at an overall growth rate of 4% annually during recent decades. A trophic mass-balance model was used to analyse the ecological carrying capacity of bivalve aquaculture in this ecoregion. The biomass of cultured bivalves is currently 13.3 t km −2 and could be increased to 62.0 t km −2 without exceeding the ecological carrying capacity. Zooplankton are a key factor limiting the ecological carrying capacity and represent a sensitive functional group within the food web system in this ecoregion. At the ecological carrying capacity of cultured bivalves in the Yellow River Estuary ecoregion, harvests would amount to 353.2 t km −2 yr −1 or a total of 4.2 million t yr −1 in this region. If the current average rate of growth in aquaculture in China is maintained, under cautious development, the biomass of cultured bivalves would reach half of the estimated ecological carrying capacity (31.0 t km −2 ) after 20 yr. This implies that there is capacity for sustainable development of bivalve aquaculture under current environmental conditions.


INTRODUCTION
Globally, China is by far the largest producer of bivalves in aquaculture. In 2018, China produced 13.4 Mt of marine bivalves of 70 species, accounting for > 80% of world production (FAO 2020). Aquaculture production increased by 4% yr −1 during the past decade. The size of the bivalve aquaculture industry can present challenges regarding coastal management, in particular in determining whether the industry is reaching or has already reached the maximum culturing capacity achievable in line with sustainable development (Costanza et al. 1997, Soto et al. 2007, Gibbs 2009). One area facing this challenge is the Yellow River Estuary ecoregion in Shandong Province, China. The region covers the Yellow River Estuary and the adjacent Laizhou Bay, and in 2018 accounted for about one-third of Shandong's 4 Mt of bivalve aquaculture production (Ministry of Agriculture and Rural Affairs 2019). Bivalve aquaculture in Shandong Province mainly comprises oyster (22%), mussel (11%), scallop (24%) and clam (33%) production.
Bivalve aquaculture management goals have advanced beyond higher production to more environmentally friendly sustainability both in China and worldwide (Gallardi 2014). Breakthroughs in artificial seeding and breeding significantly contributed to the dramatic development of bivalve aquaculture in China from the 1990s onwards (Mao et al. 2019). Longline culture and bottom culture are the primary methods of bivalve aquaculture in China; harvesting is mainly carried out manually, although mechanized harvesting is being developed (Wijsman et al. 2019).
Research into ways to deal with production breakdown caused by disease has increased (e.g. Romero et al. 2012, Travers et al. 2015, Zannella et al. 2017, and the sustainability of bivalve aquaculture has also been studied (Gibbs 2007, McKindsey 2013. Al though the culturing of bivalves influences the benthic nitrogen cycle (Kaspar et al. 1985), impacts on the benthic environment are generally low (Crawford et al. 2003, Forrest et al. 2009. Bivalves can promote benthic−pelagic coupling, making planktonic nutrients available for benthos (Smyth et al. 2013, Zhang et al. 2014, and can im prove water quality through the filter-feeding of detritus and excess nutrients and phytoplankton (Wikfors 2011, Shimoda et al. 2016. Bivalve aquaculture can also enrich biodiversity through the construction of reefs for habitat (Herbert et al. 2012, Quan et al. 2012, La Peyre et al. 2015, restoring benthic communities (Zwerschke et al. 2018) and directly providing habitats for larger macrofauna and nekton in the form of aquaculture structure (e.g. lines and cages; Tallman & Forrester 2007, Theuerkauf et al. 2022. While bivalve aquaculture has expanded, sustainability is an increasing concern. The Yellow River Estuary ecoregion is significant in terms of natural feed, spawning and habitat for various aquatic organisms (Shan et al. 2013, Meng et al. 2019, Liu et al. 2020. The Shandong Yellow River Delta National Nature Reserve is a Ramsar site (https://rsis.ramsar.org/ris/2187) located within the ecoregion, which was established in order to conserve biodiversity and ecosystem functions in the area. During recent decades, the fish assemblage structure and biodiversity in the region have greatly changed due to excessive human activity (mainly overfishing, Shan et al. 2013), and managers therefore require knowlege on to what extent the ecosystem can sustainably support aquaculture in a way that stabilizes the food web system. The issue of sustainability can be addressed through quantitatively estimating the bivalve aquaculture carrying capacity of the estuary ecoregion. The simplest way to maintain sustainability based on the metrics deemed most important to local resource managers and stakeholders is to ensure that aquaculture does not exceed its carrying capacity.
The concept of carrying capacity was derived from the logistic population growth curve, which reaches the asymptote K when the population is at maximum size (Smaal & van Duren 2019). Carrying capacity is the maximum biomass, density or number of a population that can be sustainably supported within a specified spatial area and depends on environmental factors, resources and the presence of predators, disease agents and competitors over time (Hart vigsen 2017). In bivalve aquaculture, there are 4 types of carrying capacity: physical, production, ecological and social (Inglis et al. 2000). Physical and production carrying capacity generally apply at a farm scale. Physical carrying capacity is the total farmed area that can be accommodated in a specific physical space; production carrying capacity is the stocking density that maximizes annual harvests. Ecological and social carrying capacity apply at a higher ecosystem level. Ecological carrying capacity represents the boundary conditions above which unacceptable ecological impacts occur, while social carrying capacity represents the conditions beyond which unacceptable social impacts occur.
The complex and close connection between ecological and social carrying capacity has been previously acknowledged (McKindsey et al. 2006, Gibbs 2009, McKindsey 2013, Filgueira et al. 2015, Smaal & van Duren 2019. It was concluded that ecological carrying capacity is the maximum level of change at which the impacts of bivalve aquaculture are considered socially acceptable by stakeholders (Smaal & van Duren 2019). This leads, for example, to the different ecological and social carrying capacities of developed countries and those developing countries still experiencing food security pressure (Ferreira et al. 2013). In the Yellow River Estuary ecoregion, the primary stakeholders are the local government and aquaculture farmers who wish to expand the scale of aquaculture to an ecologically acceptable level in order to decrease the amount of fishing activity while safeguarding food security. Achieving this goal requires determining the ecological carrying capacity in order to avoid disturbing ecosystem function, biomass or energy flows when planning bivalve aquaculture development in the area.
To date, only 3 previous studies have investigated carrying capacity in the Yellow River Estuary ecoregion (J. , Lin et al. 2015, 2018. J.  assessed an integrated carrying capacity index based on the pressure−state− re sponse framework for the years 2007−2009 based on indicators of coastal management effects, marine environmental quality and social support. Results revealed that carrying capacity was under heavy pressure during the period of the study, suggesting it was highly vulnerable to external impacts (e.g. increasing fishing or aquaculture activities). However, J.  did not consider ecological interaction between the indicators or identify a maximum numerical carrying capacity for sustainable aquaculture. Lin et al. (2015Lin et al. ( , 2018 evaluated the carrying capacity for shellfish and swimming crab using the mass-balance food web modelling method Ecopath (www.ecopath.org) with data collected from 2012 to 2014. However, there was only limited discussion of the ecological aspects of carrying capacity in the study area. In addition, none of this previous research specifically focussed on the carrying capacity for bivalve aquaculture. Thus, a remaining question and gap in understanding sustainable development in the Yellow River Estuary ecoregion is how much bivalve aquaculture the system can support.
The purpose of this study was, therefore, to implement a food web model for the Yellow River Estuary ecoregion using Ecopath (Christensen & Pauly 1993, Christensen 1995 to estimate the ecological carrying capacity for bivalve aquaculture. Ecopath is a massbalance, ecosystem-based food web process modelling tool that focusses on the biomass and energy flows between trophic levels. It is widely used in fisheries management for modelling a wide range of systems and management scenarios (Wolff et al. 2000, Kluger et al. 2016, Wu et al. 2016, Vasslides et al. 2017) and has been used to examine the carrying capacity of bivalve aquaculture in other systems (Jiang & Gibbs 2005, Byron et al. 2011b. As Ecopath can incorporate comprehensive trophic levels across an entire study area, including those equal to and higher than bivalves, it is more appropriate for modelling ecological carrying capacity at a regional and/or ecosystem scale than are other farm-scale modelling methods (e.g. Ferreira et al. 2007b, Newell et al. 2013. In practice, the Ecopath modelling approach relies on available biomass data for species in each functional group, as predicting biomass for unavailable functional groups via the Ecopath biomass estimation capability has been discouraged (Christensen & Walters 2004). Our analysis of the ecological carrying capacity of bivalve aquaculture within the Yellow River estuary will improve understanding of sustainable development of bivalve aquaculture in this ecologically and socio-economically important region.

Study area
The Yellow River Estuary ecoregion (12 000 km 2 ) covers the Yellow River estuary and adjacent Laizhou Bay (37°8' to 8°15' N, 118°45' to 120°15' E) and consists of a shallow-water shelf (≤15 m), with sediments consisting of soft mud and sand. The region experiences one of the strongest land−ocean interactions, producing considerable deposits and nutrient levels (Shan et al. 2013). As China's secondlongest river and the sixth longest river in the world, the Yellow River once accounted for 6% of the world's annual terrestrial sediment supply to the global ocean; however, in recent decades, it has been destabilized with dramatic reductions in water discharge and sediment load (Sun & Feng 2013, Yu et al. 2013, Shen et al. 2015. The coastal zones of the Yellow River Estuary ecoregion are economically important in Shandong Province and are experiencing rapid development of industrialization, urbanization and agri cultural practices ).

Model inputs and consideration
The food web model was constructed with 24 functional groups (Table S1 in the Supplement at www. int-res.com/articles/suppl/q014p147_supp.pdf), including 1 group comprising cultured bivalves and 3 comprising zooplankton, phytoplankton and detritus. The cultured bivalve functional group contained oysters, mussels, scallops and clams, which are the most common cultured bivalves. In total, 139 marine species were assigned to the remaining 20 functional groups based first on similarities in habitat and diet characteristics, and second, on taxonomy. All trophic levels were considered equally. Seagrass and kelp are significantly affected by the high levels of terrestrial sediment resulting in water column turbidity (Su & Tang 2002) and thus are mostly absent in this system and were not considered. Due to limits on data availability, all zooplankton species were accounted for in 1 functional group. Functional groups used are different from those in the only 2 previous Ecopath models for the study area (Lin et al. 2015(Lin et al. , 2018, possibly due to differences in data sources and research purpose. In the present study, the functional groups and their associated species were the basis for data acquisition and input parameters of the model. To model the balance between trophic levels in a food web system, Ecopath requires data for at least 3 of 4 parameters: biomass (B), production/biomass (P/B), consumption/biomass (Q/B) and ecotrophic efficiency (EE) (Christensen et al. 2005). Any one of these parameters can be calculated given in−out values for the other 3 parameters. Here values for the parameters B, P/B and Q/B were input for all functional groups. Two further necessary input components for every functional group are diet composition (DC) and fishery removals (FF, including both aquaculture harvesting and wild catch). An important but often overlooked input parameter is biomass accumulation (BA), which represents the instantaneous rate of change in biomass of a functional group. In some Ecopath models, BA has been assumed to be zero (the default value for the parameter), but this assumption leads to a dangerous bias that results in overestimates of safe harvest rates near theoretical carrying capacity (Ainsworth & Walters 2015). Detritus import was also an important parameter in this analysis, as the Yellow River brings considerable deposits and nutrients to the estuary ecoregion (Shan et al. 2013).
The Ecopath method balances the model by slightly varying the input parameters for 2 master equations (Eqs. 1 & 2) so that the mass and energy input and output are equal for each group (Christensen et al. 2005): Production = Predation + Catches + Net migration + Accumulated biomass + Other mortality (1) The variation of input parameters was restricted within the confidence limits determined through diagnostic tests (e.g. EE <1.0; production/consumption [P/C]; respiration > 0; for details, see Section 2.4). Mass and energy flows between functional groups are linked through a diet composition matrix. Some default settings were adjusted to better represent each functional group in the study area. The unassimilated consumption value was set at 0.2 for every group except zooplankton (Group 22), for which 0.3 was used to avoid underestimating the egestion of this group. All modelling was carried out using EwE 6.6.3 software (available at www.ecopath.com).

Data sources
None of the data sources used for parameterization were acquired solely for the construction of this model. The biota data acquired from external projects used multiple units of measure at different temporal and spatial scales. The various biomass units were unified by converting to dry weight standard units (t km −2 yr −1 ) and were then re-evaluated until considered appropriate for the model. The dry weight of the species sampled in situ was estimated using a conversion coefficient of 0.2 for their wet weight (Baird & Ulanowicz 1989). The basic input parameters for the model (B, P/B, Q/B, EE, BA, FF) are provided in Table 1, and DC for the groups is shown in Table 2.
B and BA for the functional groups comprising fish and other swimming species (Groups 1−8, 10 and 13−14) were estimated based on data from the years 2015−2017, acquired from a distribution and catch potential modelling study of fish in the Yellow Sea and the Bohai Sea (Zhu et al. 2020). The study collected species distribution and relative abundance data from Sea Around Us (www.seaaroundus.org) and introduced parameter values from FishBase (www.fishbase.org), SeaLifeBase (www.sealifebase. org) and the Oceanic Biodiversity Information System (http://obis.org). Based on the collected materials, Zhu et al. (2020) used a dynamic bioclimate envelope model that applied the correspondence between species' geographical distribution and habitat types (e.g. seamounts, estuaries, continental shelves) to predict current and future fish distribution and abundance throughout the Yellow Sea and the Bohai Sea, which may not consider specific ecological characteristics in the Yellow River Estuary ecoregion such as its strong land−ocean interactions.
B for the other wild species groups (Groups 9, 11−12, 15−16 and 18−23) were derived from an in situ sampling survey in 2016 (unpubl. data; Fig. 1) using otter trawls and bottom grabs. Phytoplankton (Group 23) biomass was calculated based on the chlorophyll a concentration in samples collected during the survey. BA values for the groups were estimated based on historical biomass data acquired via in situ surveys and remote sensing conducted by the National Marine Data and Information Service. Although the dataset applied different spatial and temporal resolutions to species biomass information, the Carnivorous sharks 4 Pelagic planktivorous fish 4 15 10 9 5 9 10 5 Filter-feeding sharks 6 Demersal nekton-feeding fish 8 7 7 Demersal omnivorous fish 1 8 Demersal benthos-feeding fish 1 1 9 Gobies 4  variation of the biomass during recent years (2015-2017) measured via the same methods can be used to evaluate the BA rates of the groups. B of cultured bivalves (Group 17) was estimated based on the summed production of oysters, mussels, scallops and clams in the study area in 2017 (Chen et al. 2017). The detritus biomass and annual import (Group 24) used was from the previous 2 Ecopath models in the study area (Lin et al. 2015(Lin et al. , 2018, which was based on in situ sampling surveys from 2012 to 2014. P/B, Q/B and DC were estimated based on previous studies on dynamic progress between trophic levels in the Bohai Sea and coastal areas in Shandong Province (Tang & Zhong 1990, Su & Tang 2002) and on Ecopath models of other study areas in Chinese seas (Tong et al. 2000, Cheung 2007, Cheng et al. 2009, Lin et al. 2013, 2018, Han et al. 2017, 2018. In the fish functional groups (Groups 1−9), P/B and Q/B were updated using the life history tool in FishBase (www.fishbase.org/, accessed 13 April 2020), based on the physiological and environmental features (e.g. fish length and sea temperature) of the samples collected during the in situ survey con ducted in 2016 (Fig. 1). Data for crustacean functional groups (Groups 11−14) Q/B was updated using the species dataset Sea Around Us (www.seaaroundus.org/, accessed 13 April 2020). Studies on the diet composition of marine species in Chinese seas (Yang 2001a,b, Xue 2005, B. Zhang et al. 2007, 2012, Zhang 2018) were used to optimise the estimation of DC. As the fish assemblage structure in the study area has changed greatly during recent de cades (Shan et al. 2013), newer published studies were prioritised in the estimation to better represent current diet composition. The prey of cultured bivalves (Group 17) was set to be the same as that for wild bivalves (Group 16), since local wild bivalve species were originally chosen to be cultured in this estuary ecoregion and included both filter feeders and deposit feeders (Ministry of Agriculture and Rural Affairs 2017). Cultured bivalves (Group 17) were set to be not preyed on by any other groups (i.e. no predators), following previous research using the Ecopath approach for estimating carrying capacity of bivalves/shellfish in other study areas (e.g. Jiang & Gibbs 2005, Byron et al. 2011b). The setting is used to conservatively evaluate carrying capacity and maximize potential impacts from expanding bivalve aquaculture. The premise of minimal predation on cultured bivalves may also replicate conditions of ideally managed aquaculture farms and yet allow these bivalves to prey on lower trophic-level groups.
Data for removal by the fishery fleet (aquaculture harvest and/or wild catch) within the study area in 2016 originated from the China Fishery Statistical Yearbook 2017 (Ministry of Agriculture and Rural Affairs 2017) and the Shandong Statistical Yearbook 2017 (Chen et al. 2017).

Diagnostics
Before balancing the model, the validity of the input parameters was evaluated using the diagnostic test method applied by Link (2010) and Byron et al. (2011b). The diagnostic testing was independent of the Ecopath algorithm assumptions and was conducted prior to the Ecopath automated balancing routine, providing a manual check of the validity of the data sources and greater control over the massbalancing. Diagnostic testing evaluates the cohesiveness of variously sourced data under multiple spatial and temporal scales, revealing any discrepancies, allowing the model to be holistic rather than piecemeal (Heymans et al. 2016). Diagnostic tests were performed on each functional group in sequence of decreasing trophic level; the factors tested were biomass, the ratio of biomass to total primary production, vital rates (production [P], consumption [C], respiration [R]), ratios of vital rates (P/C, P/R) and total consumptive removals (Fig. 2). The following validation rules (Byron et al. 2011b) were applied: (1) biomass and vital rates should decrease with trophic level; (2) ratios of biomass to primary production should be <1; (3) vital rate ratios should be within ecologically acceptable limits (i.e. P/C < 0.5); (4) consumption should be higher than production within each functional group; and (5) total consumptive removals should be lower than production. When the diagnostic tests indicated a failure to comply with a rule, input parameters were adjusted for ecological integrity and validity by a multiple of 2 / 3 , 1 / 2 , 3 / 2 , 2, 3, 10 or 25 as necessary.
Further improvements of input parameters were made and data uncertainty was addressed using the 153 Fig. 2. Pre-balancing diagnostics on functional groups ordered by trophic level. See Table 1 for group numbers. (a) Trophic decline in biomass required across species groups. (b) Ratios of biomass (B) and production (P) to total primary production (PP), all required to be below 1. (c) Vital rates required to decline with trophic level, with higher consumption than production and respiration. (d) Vital rate ratios: consumption/biomass (Q/B), respiration/biomass (R/B) and production/biomass (P/B). (e) Production/consumption (P/C) required to be lower than production/respiration (P/R). (f) Total consumptive removals required to be less than production and consumption Ecopath tool 'Pedigree'. This tool evaluates statistical uncertainty by giving each parameter a confidence interval and ranking these intervals by their respective data sources (e.g. data sourced from sampling is ranked higher than a modelling estimation). A probability distribution for each parameter was then established, based on the confidence intervals, through the Monte-Carlo parametrization carried out by the 'Ecoranger' module (Christensen et al. 2005).

Model outputs
To evaluate influences between the functional groups in the food web system in the study area, mixed trophic impact analysis in Ecopath was used to determine how much a group could be affected by slight perturbation in another group, based on their food web characteristics. The analysis involves calculating the impact of an infinitesimally small increase in the biomass of one group on other groups (Christensen et al. 2005), which can be output by Ecopath as a matrix of impacts between groups (i.e. a mixed trophic impact plot).
Ecopath was also used to output the balanced model statistics that characterized the system, providing informative system measures such as throughput, cycling index and pathways. Total system throughput represents the scale of the system in terms of mass and energy flow (Ulanowicz 2012). It sums all flows in a system, including consumption, ex port, respiration and flows to detritus (Christensen et al. 2005). The cycling index is the fraction of energy throughput in a food web system that is recycled (Finn 1976), which correlates with system maturity, resilience and stability (Christensen et al. 2005). The number of pathways reflects the redundancy and stability of the system.

Assessment of carrying capacity
Ecological carrying capacity was analysed following the methods of Jiang & Gibbs (2005) and Byron et al. (2011b). Starting from the balanced model, the biomass of cultured bivalves (Group 17) and those harvested (i.e. their landings) were increased while keeping all other parameters constant until the model became unbalanced and no longer represented the current condition of the system. The level of biomass just prior to any change in the model was the ecological carrying capacity of cultured bivalves.
When the model was balanced at the ecological carrying capacity of cultured bivalves, the robustness of the biomass parameter was examined by changing an individual biomass value from each group by factors of 0.01, 0.5, 2, 10 and 100 of their original value, while keeping all other values constant. The robustness of the functional group's biomass with respect to perturbations at carrying capacity was indicated by how much the biomass value could vary while the model remained balanced.
Production carrying capacity analysis carried out by removing specific functional group(s) (e.g. Byron et al. 2011a,b) was not performed in this study, as the removal of any functional group(s) (such as a group disappearing because of inadequate food) would cause a fundamental change in the ecological structure in the study area, including a potentially significant decline or increase in primary production, and leading to risks in food security, which were not acceptable to stakeholders.
The Pedigree index score was 0.758 with a measure of fit of 5.319 (Table S2), demonstrating high data source quality (Morissette et al. 2006). Confidence was high with respect to biomass input parameters and low with respect to catch (landings), while the input parameters for sharks, crabs, polychaetes and detritus were lower than in the other groups, as their biomass differed greatly from the expected trendline (Fig. 2a).
Mixed trophic impact measures the direct or indirect impact of a functional group on another group, based on food web characteristics (Fig. 4). Detritus has a large positive impact on the pelagic omnivorous fish and zoobenthos groups. The largest negative impacts were that of Scomberomorus niphonius on demersal omnivorous fish and on itself. The shrimp group generally had negative impacts on zoobenthos and on itself. Zooplankton also had negative impacts on lower trophic level groups. Changes in biomass of the cultured bivalve group evenly affected every other group, including negative feedback on itself.
The total number of links between functional groups was 33 709, with a mean pathway length of 2.8. There were also 154 pathways that started from a group and returned to it. The cycling index (Finn 1976) was 4.5%, with a mean path length of 2.5. The total throughput (i.e. the total energy in the food web system) was 6210 t km −2 yr −1 (Table S3). Only 0.35% of the total throughput, excluding detritus (i.e. predatory cycling index), was cycled (7.8 t km −2 yr −1 ). The largest energy flow was through primary producers (Group 23), which accounted for 39% of the total throughput (Table S3). Zooplankton consumed 77.7% of the primary production and were simultane-  (Table 2) and in direct competition with the cultured bivalves for food (Fig. 4). Most of the energy in the food web was cycled through producers or diverted to detritus (32%).

Carrying capacity
The biomass of cultured bivalves was 13.3 t km −2 and could be increased to 62.0 t km −2 without exceeding the ecological carrying capacity (Table 3). Under the current average rate of growth of bivalve aquaculture in China (4% yr −1 ), bivalve biomass would reach half of the estimated ecological carrying capacity after 20 yr. Given the total area of the Yellow River Estuary ecoregion, the total potential ecological carrying capacity for cultured bivalves is 0.74 Mt. At ecological carrying capacity, the landings (harvest) of cultured bivalves could reach 353.2 t km −2 yr −1 , or 4.2 Mt yr −1 . The harvest could thus be more than 6 times greater than that currently reported in the study area.
In the estimation of carrying capacity, zooplankton production was totally consumed (EE ≥ 1), causing the model to become unbalanced when the biomass of cultured bivalves exceeded ecological carrying capacity in the study area (Table 3). Zooplankton were thus the limiting factor when the Yellow River Estuary ecoregion was at ecological carrying capacity for cultured bivalves. In the robustness analysis (Table 4), a slight change in zooplankton biomass led to an unbalanced model and unacceptable change in the food web. At ecological carrying capacity, the biomass parameters of swimming species were more robust in the model, especially those species with very low biomass (sharks). The model also remained balanced when the biomass of Oratosquilla oratoria, a popular fished species in this ecoregion, was halved. Phytoplankton and detritus were similarly highly robust at the ecological carrying capacity of cultured bivalves.

DISCUSSION
Cultured bivalves are an important part of the food provision system in China and in the Yellow River estuary study area, with a rapid increase in the bivalve aquaculture industry during recent decades (Ministry of Agriculture and Rural Affairs 2019). Compared with previous Ecopath modelling that did not include cultured bivalves (Lin et al. 2015(Lin et al. , 2018, this study documented clear differences in diagnostics, automated mass-balance routine and final model outputs, and in particular showed that energy cycling and system activity in the model that included cultured bivalves were much higher than in those that did not. In the current model, the cultured bivalve functional group contributed 62% of throughput excluding plankton and detritus (Table S3). Impact analysis showed that cultured bivalves evenly affected almost all of the other functional groups, including themselves (Fig. 4). Unlike wild bivalves, cultured bi valves were not allowed to be preyed on by other groups nor to prey on other groups except for plankton and detritus. This resulted in cultured bi valves being the primary competitor with other functional groups using primary production and low-trophiclevel groups as a source of food. The cycling index (Finn 1976) of this model (4.5%) is similar to those of the previous 2 models. The cycling index was 3−6% in the model based on data acquired during 2012−2013 (Lin et al. 2015), and 2.8% in the model based on data acquired during 2013−2014 (Lin et al. 2018). The low cycling index, which represents a fraction of energy retained (Christensen et al. 2005 Table 3. Changes in the Yellow River Estuary ecoregion Ecopath model while estimating carrying capacity using the automated mass-balancing routine. Numbers in bold type are the calculated ecological carrying capacity biomass values for cultured bivalves (Group 17) and their landings. Landings are the predicted amount of harvest when the ecotrophic efficiency (EE) of cultured bivalves equals 0.95. Any biomass value below the carrying capacity did not affect the balanced model. Any biomass value above the carrying capacity resulted in an unbalanced model as it led to the EE of another group or groups becoming ≥ 1, signifying that the group(s) were completely consumed length in our model (2.5), suggests that energy is not efficiently recycled and materials are not retained over time in the Yellow River Estuary ecoregion. In addition, as the total throughput represents the size of the food web's energy flow (Ulanowicz 2012), the low total throughput of our model (6210 t km −2 yr −1 ) indicates a low level of energy flow through the food web system in the Yellow River Estuary ecoregion. The low level of throughput may be due to low levels of photosynthesis. The physiology of benthic plants (e.g. seagrass, kelp) in the study area was greatly affected by high levels of turbidity resulting from the large amounts of terrestrial sediment brought to the region by the Yellow River (Su & Tang 2002). The throughput value was about double that observed in 2 previous Ecopath models (Lin et al. 2015(Lin et al. , 2018, in which the biomass of functional groups was generally lower than in this study (e.g. the biomass of zooplankton was estimated as 2~5 t km −2 whilst it was 10.9 t km −2 in the current study). Increasing the number of bi valves cultured could improve energy cycling and system activity by strengthening benthic− pelagic coupling (Stolpovsky et al. 2018, Ricci et al. 2019), because bivalves filter materials from waters and then excrete organic wastes that are physiologically available to benthic consumers and plants.
The estimated ecological carrying capacity of cultured bivalves within the Yellow River Estuary ecoregion (62.0 t km −2 ) is lower than that in other estuary and/or bay regions, despite the considerable deposits and nutrients brought to the area by the Yellow River every year (Shan et al. 2013). It is similar to the carrying capacity estimated in the model for the Tasman and Golden Bays in New Zealand (65 t km −2 , Jiang & Gibbs 2005), but lower than that for Narragansett Bay in the USA (297 t km −2 , Byron et al. 2011b). It is likewise lower than the input biomass parameters of wild bivalves used in Ecopath models for some other study areas (e.g. Arbach Leloup et al. 2008, Zajac et al. 2008. The moderate value is primarily due to an insufficient biomass of zooplankton (Group 22) as a food source for cultured bivalves (Table 3, Fig. 3), and partly due to zooplankton competing with cultured bivalves for consumption of phytoplankton and detritus. Although all zooplankton species were accounted for as 1 functional group, bivalves usually  (Forrest et al. 2009) and do not directly compete with other large zooplankton consumers (pelagic fish, filter feeding sharks, etc.) which consume meso-zooplankton. Therefore, the biomass of zooplankton actually available as food for cultured bivalves in the study area could be less than that represented in the balanced model output, which would also imply a lower carrying capacity in practice. The biomass and energy flow through other functional groups were too small (Tables 1, Fig. 3) to affect the cultured bivalves and their carrying capacity (Fig. 4). As a food source, detritus may be another significant factor limiting carrying capacity for cultured bivalves in the study area (Table 3, Fig. 3).
Researchers have estimated that cultured shellfish can consume 4 times as much detritus as phytoplankton (Ferreira et al. 2007a). The low level of photosynthesis in the food web system might also generally limit the carrying capacity in the study area.
As the current biomass of cultured bivalves (13.3 t km −2 ) is less than a quarter of the estimated ecological carrying capacity (62.0 t km −2 ), there is no indication that the bivalve aquaculture industry is currently food-limited, nor that bivalves are altering the food web structure in the Yellow River Estuary ecoregion at the current average rate of growth of this industry (4% yr −1 ). Bivalve aquaculture in the study area continues to increase (Guo et al. 2019), but potential production biomass (landing) output from the Ecopath model (353.2 t km −2 yr −1 ) is much higher than current production (54.3 t km −2 yr −1 ) (Table 3), so bivalves are not being farmed at maximum potential. This could be due to underdeveloped artificial seeding and less mechanized harvesting in past decades (Mao et al. 2019, Wijsman et al. 2019. The robustness test (Table 4) shows that when the food web system approaches the ecological carrying capacity for cultured bivalves, a slight increase beyond carrying capacity would have significant ecological implications. The biomass of sharks and fish would be better tolerated as the biomass of cultured bivalves advances towards its carrying capacity (Table 4). Although inadequate levels of zooplankton limit the carrying capacity of cultured bivalves, increasing zooplankton biomass results in an unbalanced model. Due to competition between zooplankton and cultured bivalves for prey, an increase in zooplankton could result in inadequate detritus and/or phytoplankton so that the model fails to balance.
Zooplankton may play an important and distinct role in the food web system of the Yellow River Estuary ecoregion, as in most other estuary ecoregions the ecological carrying capacity of cultured bivalves is not limited by zooplankton biomass. However, simply increasing zooplankton biomass by any means (which may not be technologically feasible) could cause a significant change in the food web system as it would lead to detritus and/or phytoplankton being completely consumed (Table 4). While the system is sensitive to zooplankton biomass, the biomass accumulation (BA) rate of zooplankton in the model (+ 0.22) indicates a large seasonal change in zooplankton biomass (Table 1). In addition, the cephalopod and echinoderm functional groups had low robustness and high BA rates, while phytoplankton had a BA rate of −0.56 and its robustness was low when biomass was decreasing (Table 4). These factors contribute to the vulnerability of the food web in the study area. The new trend of integrated aquaculture ('eco-farming') in China (Wijsman et al. 2019), such as shellfish−fish, shellfish−shrimp and shellfish−seaweed, might help to strengthen the general robustness of the entire food web system by better accounting for trophic interactions. This type of integrated multi-trophic aquaculture has been conducted in Sanggou Bay, China, for several years, with over 30 species including bivalves producing a total of > 240 000 t of seafood each year from approximately 100 km 2 of production space (Fang et al. 2016).
The model outputs suggest that bivalve aquaculture does not significantly alter food web function and is at or below its ecological carrying capacity, at least as modelled here at current biomass of cultured bivalves (Table 4). This information can be used to develop aquaculture management in the Yellow River Estuary ecoregion by balancing ecological conservation and socio-economic considerations. In addition to the stakeholders' primary concerns for safeguarding food security, growth of the bivalve aquaculture sector could also provide other services not examined here like providing habitat for other organisms (Zwerschke et al. 2018), mitigating eutrophication and enhancing nutrient cycling via shell and tissue production (Carmichael et al. 2012). A new data series for the biomass of each functional group (including some groups lacking in this model like separate zooplankton categories, and bacteria which could also be preyed on by bivalves and sensitive to environmental conditions) could be collected over the next few years, and analysed via the Ecopath method (or Ecosim if the data series are temporally continuous and complete) to estimate potential temporal variations in the ecological carrying capacity in the study area. Such an expanded modelling effort could simulate possible environmental changes in terms of further nutrient inputs derived from development on land, and/or changes in wild species' distributions caused by global climate change. Information for each functional group at a seasonallevel temporal resolution could be introduced in the estimation of ecological carrying capacity for a better representation of the ecological conditions in a typical temperate estuary region.
It should be noted that the modelled ecological carrying capacity is a theoretical limit only, is based on the accurate assessment of trophic contributions made by separate functional groups, and should be treated with caution in practice. In many fisheries management systems, half the carrying capacity is an acceptable target for maximum sustainable production (Mace 2001). With a target of half the ecological carrying capacity of cultured bivalves (31.0 t km −2 ) and assuming an industry growth rate equal to China's current national average of 4% yr −1 , it would take 20 yr for bivalve aquaculture in the study area to approach maximum sustainable production. This implies that changes to the present plan for development of the bivalve aquaculture industry are not necessary. This target of achieving half carrying capacity provides sufficient operating space to address differences between the theoretical model and practice, and any other variables that could potentially affect carrying capacities, such as factors affecting physical carrying capacity (e.g. hydrodynamics, water quality and available space), production carrying capacity (e.g. stocks of wild or cultured competing species) and social carrying capacity (e.g. food security, market; Smaal & van Duren 2019).