Implementing a finite-volume coupled physical-biogeochemical model to the coastal East China Sea

Several models for estuarine physical processes and biogeochemistry have been developed over last decades. One of the most comprehensive coupled model systems, Finite Volume Community Coastal Model 10 (FVCOM) coupled with European Regional Seas Ecosystem Model (ERSEM) through the Framework for Aquatic Biogeochemical Models (FABM) has been implemented to a high resolution coastal East China Sea (ECS), which encompassed complex coastal zone and part of continental shelf. Physical model was assessed by traditional univariate comparisons, while a rigorous model skill assessment was conducted for coupled biological model. The model system’s ability to reproduce major characteristics both in physical and biological environments was 15 evaluated. The roles of physical, chemical and environmental parameters on the biogeochemistry of the ECS were extensively studied. This work could form a significant basis for future work, e.g. the response of biogeochemical flux to physical mechanism.


Introduction
Extensive ocean models have been developed over last several decades to serve as tools for research and maritime projects. A demand for explicit modeling of combined physical, chemical and biological systems begins on a growing realization that biogeochemical state cannot be inferred from their physical properties alone (Blackford et al., 2004;Tomasz et al., 2014). Traditional numerical models have been 45 constructed based on simplified assumptions on the functionality of complex marine ecosystem. Most of them failed to simulate important biogeochemical processes, because the models did not consider essential features, such as explicit carbon cycling, microbial food dynamics, the role of key functional groups and multiple nutrient limitation to primary production (Mateus et al., 2012). In recent years, marine ecosystem models have been explored to understand, quantify and estimate biogeochemical 50 processes in seas and oceans. These models vary in complexity from simple four-compartment Nitrate, Phytoplankton, Zooplankton, Detritus (NPZD) pelagic models (Oschlies et al., 2000;Dabrowski et al., 2013) to more complex multi-functional group models describing ocean biogeochemistry and lower trophic food web (Mateus, 2012;Flynn, 2010;Follows et al., 2007;Wild-Allen et al., 2010). European Regional Seas Ecosystem Model (ERSEM) is one of the most established and complex ecosystem 55 models for lower trophic levels of marine food web in use, which assesses over 40 state variables with benthic-pelagic process (Baretta et al., 1995). The model has been applied in a wide number of contexts that included short-term forecasting (Edwards et al., 2012), ocean acidification (Blackford and Gilbert, 2007), coupled climate-acidification projections (Polimene et al., 2014), and biogeochemical cycling (Wakelin et al., 2012). 60 the economy. Previous research have focused on this ecosystem hazards by illustrating lower nutrient 70 food web along coastal region of ECS (Wang et al., 2006;Guo et al., 2014;Li et al., 2008;Ye et al., 2015;Zhang et al., 2004). However, most of these studies were based on either limited field surveys or simplified laboratory experiments. There were also defects in providing environmental drivers offline, the simple zero-dimensional box biological model, and limited spatial and temporal coverage (Wang et al., 2013;Zhu et al., 2005). Online three-dimensional physical and biogeochemical fluxes need to be 75 considered for more realistic representations.
In this study, we implemented the high-resolution FVCOM coupled to the ERSEM through the FABM framework. The paper was organized as follows. Section 2 briefly described the study area and observations used. Section 3 focused on the FVCOM-ERSEM model, and specific setup. Model skill assessment was carried out in Section 4. A discussion was further explored in Section 5. Conclusions 80 were drawn and future work was discussed in Section 6.

Study area and observations
Study area covers the coast of ECS extending from southern Taiwan Strait to northern Yangtze River (YR) water system. Model domain is 117 -124.5°E and 22-33°N within a 174 m-isobath (Fig.1).
Physical environment of ECS has a distinct seasonality feature at mid-latitude and influenced by 85 anthropogenic stresses from adjacent landmass, as well as mixing from several principal water types.
From the north to the south, the domain comprises several important estuaries, including: YR, QR river, OR, Aojiang River (AR), Feiyunjiang River (FR), MR and Jiulongjiang River (JR) and empties into ECS. The topography of model domain is derived from ETOPO1 (1-minute gridded data) for the open ocean region (Amante and Eakins, 2009) and nautical chart for estuarine areas. 90 Monitored data were collected from 5 ecological buoys monitoring at the Zhejiang coast in May 2019, including temperature, salinity, conductivity, pH, dissolved oxygen (DO), dissolved oxygen saturation, turbidity, chlorophyll-a (Chla), phosphate, nitrate and ammonia. Meanwhile, tidal level observations were also collected at several tidal stations to evaluate model performance.

Hydrodynamic model
The numerical model used in this study is unstructured grid based, free-surface, 3-D primitive equations Finite-Volume Community Ocean Model (FVCOM) ocean model described in detail by Chen et al. (2003a). To date, current version is fully coupled ice-ocean-wave-sediment-ecosystem model system 105 with options of various turbulent mixing schemes, generalized terrain-following coordinates and wet/dry treatments. Finite-volume approach combines finite-element method for geometric flexibility and finitedifference method for simple discrete structures, in order to enhance the computational efficiency.
Multiple dynamical forces, including river runoff, astronomical tide, waves, mean flow, wind, etc. and seasonal temperature, salinity and density, coexist and interact in study area. Therefore, unstructured, 110 finite-volume ocean model can fit to ECS situation sensationally.

Biogeochemical model
ERSEM was developed in the 1990's, which is one of the most established ecosystem models for the lower trophic levels of the marine food web (Butenschön et al., 2016). The current model release https://doi.org/10.5194/os-2020-47 Preprint. Discussion started: 9 June 2020 c Author(s) 2020. CC BY 4.0 License. contains the essential elements for the pelagic and benthic parts of the marine ecosystem, including the 115 microbial food web, the carbonate system, and calcification.
Trophic structure is defined on the basis of a predatory action of consumers on producers, bacteria, and themselves. A benthic module is implemented to estimate the mineralization of sinking organic matter, nutrients and oxygen flux at bed-water interface. A pelagic model comprises more than twenty-two major state variables: light, producers, consumers, decomposers, pelagic organic matter (dissolved labile, 120 semi-labile and semi-refractory DOM, small-size, medium-size and large-size POM), benthic organic matter (dissolved, particulate and refractory), nutrients (nitrate, ammonium, phosphate and silicate), and oxygen. For this simulation, we have considered four types of primary producers: diatoms, nanophytoplankton, picophytoplankton and microphytoplankton, and three group of consumers: mesozooplankton, microzooplankton, and Heterotrophic flagellates. 125 The FABM enables complex biogeochemical models for marine and freshwater systems to be developed as sets of stand-alone or process specific modules (Bruggeman & Bolding, 2014). It has been coupled to many hydrodynamic models including GOTM (General Ocean Turbulence Model), ROMS (Regional Parameter values for each generic type model were listed in Tables 2 to 4. Whenever possible they have been adopted from original study performed in similar complexity (Butenschön et al., 2016). Some parameters like maximum specific productivity, mimimal specific lysis rate, assimilation efficiency of mesozooplankton were estimated from references of ECS (Guo et al., 2014;Li et al., 2008;Gin et al., 135 1998;Wang et al., 2006).

Model configuration
The computational domain was divided into 102, 688 non-topped triangular cells with 53, 512 grid nodes. The resolution at open boundary was set up to 15 km, and refined to ~200 m around riverine channel. 140 The model was forced by realistic tide, river discharge and atmospheric conditions. The tidal forcing was imposed using the TPXO7.2 data (Egbert and Erofeeva, 2002), which provides 8 primary harmonic constitute to predict ocean tide. Inputs for fresh water were prescribed using climatological monthly https://doi.org/10.5194/os-2020-47 Preprint. Discussion started: 9 June 2020 c Author(s) 2020. CC BY 4.0 License. values for respective main channels (Fig. 2), derived from China Water & Power Press (http://www.waterpub.com.cn/). Riverine inputs of salinity were set to 0 psu. The temperature and 145 suspended sediment concentration were collected according to multi-year averaged monthly data sets (Editorial Board of China Bay Survey). The surface wind forcing, heat flux and precipitation/evaporation were acquired from 6-hour Reanalysis data of NOAA/s National Centers for Environmental Prediction (NCEP) (ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis2/gaussian_grid).
Initial conditions for the temperature and salinity were derived from GDEM (Generalized Digital 150 Environmental Model). Open boundary conditions for temperature and salinity were extracted from an ocean reanalysis data set of SODA (Simple Ocean Data Assimilation).
Surface boundary conditions were prescribed as no-flux for all biogeochemical state variables. Monthlyaveraged nutrients were imposed at riverine boundary of YR in Table 1, including nitrate nitrogen, ammonia nitrogen, phosphorus, silicate (Wang, et al., 2013;Liu, 2002;Xu, 2019). Yearly-averaged 155 nutrients were specified at other riverine boundaries. The Nitrate nitrogen, ammonia nitrogen, phosphorus, silicate concentrations were set as 80.4, 2.26, 1.53 and 120 μmol l-1 respectively at QR; the corresponding concentrations were set as 20, 5.5, 1 and 150 μmol l -1 at OR, AR and FY; and 53.3, 13, 0.5 and 221 μmol l -1 at JR. Initial and open boundary conditions for the ecological model properties (phosphate, nitrate, oxygen and silicate) were derived from WOA09 (World Ocean Atlas 2009), and the 160 Chla was from OC_CCI (Ocean Colour Climate Change Initiative). The initial values of ammonium were given as a homogeneous constant value, 1.0 mmol N m -3 .

165
An initial period of 1 month, January 2018, was used as a spin-off period for tide in barotropic mode, and the model was run for 11 months until December in baroclinic mode with an initial tide fields. The

Model skill assessment
This section concerned physical and biological response of the model during numerical period.
Traditional univariate comparisons were used to assess physical model skill. Time-series comparisons 180 of water level, sea surface temperature, and sea surface salinity were presented in Figs. 3-5, whereas Table 5 showed the comparisons of harmonic constants for 8 main astronomical constituents.
A rigorous model skill assessment was conducted for coupled biological model, thus the model's capabilities were tested and understood. Herein, additional approaches were explored to validate the complex biological model performance, e.g. the Percentage of Bias (PB) and the Adjusted Relative 185 Mean Absolute Error (ARMAE). We also attempted multivariate comparison of the modeled and the observed using the Cost Function (CF) to minimize model-data misfit.

Physical fields
Graphically comparing the modeled with the observed could be a useful way to assess model performance. We plotted time-series water level, SST and surface salinity. To quantify the differences 190 between data and modeled results, the root-mean-square (RMS) error ( ) and the correlation coefficient ( 2 ) were employed as major skill assessment index of physical fields to indicate average discrepancy (Taylor, 2000).
where and were for the modeled and the observed, respectively. ̅ and ̅ were mean values, and and were standard deviations of and , respectively.

Tidal analysis
In coastal and estuarine environment, tidal current represented one of main forcing on biogeochemical dynamics. It was crucial to correctly simulate tide propagation along the coast that could be 200 characterized by water surface elevation. We selected six tidal stations to validate temporal water level of May 2019 (Fig. 3). As shown in Fig.3, simulated water levels were in good agreement with observations over entire measurement period except for the RA and the PT. The correlation coefficients 2 ranged from 0.94 to 0.99, and the root-mean-square (RMS) error from 0.155 to 0.59 m. The R station lay near AR estuary, where the coastline and terrain have been modified due to human 205 reclamation in past several years, thus it appeared that flood tide could not pump high enough. The PT located near eastern coast of the PT Island, where the resolution of the terrain was quite coarse.
Harmonic analysis were also carried out using the T_Tide Toolkit with observed tidal level and computed water level. Harmonic constants of eight main tidal constituents were listed in Table 5.
Generally, maximum deviation of the amplitude was less than ~10 cm, and phase difference was less 210 than ~20 degree. Compared to tidal ranges of coastal areas, model and observational data are in good agreement with each other.

Surface sea salinity and temperature
Temperature structure was one of major limiting factor on primary productivity in coastal ECS. The  (He et al., 2007;Chin et al., 1965). Salinity influenced penetration pressure of alga, thus physiological state would be changed to some extent. Each group of phytoplankton had suitable specific range of the salinity, e.g. prorocentrum donghaiense Lu prefers 25-35 psu (He et al., 2007). Therefore, 3-D baroclinic fields were significantly important for the modeling for biochemical cycling.

245
Compared to physical oceanography, field observation for chemical and biological oceanography was scarce and this remained an obstacle to improve coupled biogeochemical model system (Ward et al., 2010). The Chla was commonly collected in the estimation of coupled biogeochemical model for wide availability of the data, both in-situ and remotely sensed, and was a focus of this model assessment.
Validation for nutrients and Chla at observed sites were also collected. offshore region, about twice than simulated values (Fig.8). The simulated phosphate was at good 255 agreement with the observations (Fig. 6). Generally, a reasonable agreement for nutrients was achieved in terms of concentration magnitude. The wz01, located near the coastline of the Leqing Bay, was affected by pollutants input. The model underestimated the observed NO3 at this site, which approximately equaled to the initial fields.

Chla and DO
Figures 9 and 10 presented the simulated Chla and DO at four stations (nb03, zs03, wz01 and wz02).
The range of oxygen concentrations at each station were well reproduced, with some overestimation at the wz01. Overall results suggested that oxygen budget in the system was satisfactorily achieved. 270 A comparison with Chla showed that the model has broadly captured observed variation, although with lower magnitude of the nb03 and wz02, where bloom peak occurred in late-spring month. Field values could reach to 15.0 mg Chla m -3 at nb03, and almost 35.0 mg Chla m -3 at wz02. Figure 11 showed the comparison between modeled Chla and ocean color product of the Sea WiFS. Simulated Chla were in https://doi.org/10.5194/os-2020-47 Preprint. Discussion started: 9 June 2020 c Author(s) 2020. CC BY 4.0 License. typical average range of 0.0-20 mg Chla m -3 with an overestimation near the YR offshore areas, and an 275 underestimation along the Zhejiang Province coast, e.g. OR estuarine areas. Generally, the model was able to reproduce major temporal and spatial variability, although there was a mismatch existed in late spring abundance. The reasons were possibly due to followings: (1) input values of properties of coastal rivers used were not daily but monthly averaged, (2) inaccurate parameters for chlorophyll synthesis, growth rates, etc., (3) absent of suspension of sediments and frequent high 280 turbidity in tidal estuaries along the coast.

Quantitative assessments
Following statistical measures: the Cost Function (CF), the Percentage of Bias (PB) and the Adjusted Relative Mean Absolute Error (ARMAE), were introduced to assess biological models usefulness as 290 tools in decision-making process. These statistics delivered model performance were defined as follows (Tomasz et al., 2014):  305   Table 6. The model was assessed by CF as excellent/very good for ammonium, phosphate and Chla, and good for nitrate and DO. The values of ARMAE represented the relative error over and above the estimated error in the observations. The model scores good for ammonium and Chla, excellent/very good for phosphate and DO, and poor for nitrate. The assessment based on PB was more rigorous, as PB showed the bias normalized by the observations rather than standard deviation. The model was 310 classified as excellent/very good for DO, and poor/bad for the remaining state variables.  Figure 12 showed modeled distribution of the nutrients in May of bloom peak month. Nutrients exerted the control to some extent on phytoplankton composition in some systems. It was particular relevant for silica because the model estimated a decrease of this nutrient from the inshore region to offshore region.
Both nutrients showed a higher value near estuarine and nearshore areas, especially the YR estuary and 320 coastal areas of Zhejiang Province. Affected by the interaction of the YR diluted water and longshore current, four type of nutrients roughly kept similar trend that decreased from the inshore plume region to offshore region. Nutrient front was formed offshore with different patterns. Nutrients were depleted in the well-mixed areas away from the front. The patterns appeared roughly same trend with observations in May 2015 (Ye et al., 2015), although the data were collected in a different year. 325

Phytoplankton
Current pelagic model of ERSEM comprised four functional types for primary producers based on single trait size, with the exception of the requirement of the silicate by diatoms and an implicit calcification 330 potential of nanoflagellates. This lead to the four classes of diatoms, and pico-, nano-, and microphytoplankton (Butenschön et al., 2016). As shown in Fig. 13a, high values of diatoms appeared at the QT and OR river mouth, and mean value of 2-5 mg C m -3 in plume region and offshore areas. The nanophytoplankton occupied almost entire domain except for upstream river mouth with relatively smaller values (Fig. 13b). For picophytoplankton, obvious higher values occured outside the YR mouth, 335 where Guo et al. (2014) also concluded the similar pattern according to two curises in August 2009 (Fig.   13c).

Zooplankton
Three predators were considered including the microzooplankton, mesozooplankton and heterotrophic flagellates according to their size classes. They predated on different prey types, including cannibalism (Butenschön et al., 2016). Grazing was treated as major control on phytoplankton abundance in pelagic system, but the lack of field data on trophic structure in coastal areas limited assessing the magnitude of 345 this control in coastal ECS. Observational data were not available for the model, moreover, zooplankton research was outside the aim of our present study, so we only plotted spatial distributions of the three zooplankton for the integrity and made some rough estimates (Fig. 14). It showed clearly that heterotrophic flagellates occurred mainly outside YR estuary and along the coast of ZP. High values of microplankton appeared near offshore plume areas, while the mesozooplankton grew along the coast 350 and high value could be seen in offshore plume region. Apart from these runoff, there are many small rivers or streams, which affected circulation, especially https://doi.org/10.5194/os-2020-47 Preprint. Discussion started: 9 June 2020 c Author(s) 2020. CC BY 4.0 License.
water quality near river mouths. Unfortunately, most of those small rivers are not continuously measured. 360 Biological state variables of rivers including nutrients, DO and biomass, are much more difficult to be collected, which usually are assumed with historical data or references. Accurate field data of river mouths would improve the accuracy of model simulation greatly, particularly around estuarine areas.
Model results are dependent on the parameterizations of physical and biological processes in component modules of FVCOM-ERSEM. For example, number and type of nutrients influenced biological results. 365 Adsorption of suspended sediment and resuspension near bed affect the nutrient cycle. Also the dynamics of algal groups determined lower trophic food web.

Conclusions and future work
This paper presented a 3-D finite-volume physical-biogeochemical coupled model of coastal ECS. We implemented the 3-D baroclinic physical model FVCOM, which utilized a triangular horizontal grid to better fit estuarine and coastal geometry, coupling with ERSEM, a well-established ecosystem model 375 for lower trophic levels of marine food web, through the framework of FABM.
The model performance was assessed by extensive validation for major characteristics both in physical and biological environments, including the variables of water elevation, temperature, salinity, surface Chla and nutrients concentrations. Due to the limitation of observational data, we evaluated simulated results of phytoplankton on a qualitative basis. The model was capable to reproduce main observed 380 temporal and spatial features for phytoplankton and nutrients. The nutrients and phytoplankton distributions of coastal region of the ECS were discussed briefly in this study. This integrated ecosystem model could be further explored to assess the response of biogeochemical process to physical, chemical and environmental parameters for coastal ECS. Certainly, we will continue to improve the physical model and biogeochemistry model parameter space. 385