The second-generation real-time ecological environment prediction system for the Guangdong–Hong Kong–Marco Greater Bay Area: Model setup, validation, improvements, and online visualization

With the rapidly growing population and socioeconomic development of the Guangdong–Hong Kong–Marco Greater Bay Area of China, inputs of diverse contaminants have rapidly increased. This poses threats to the water quality of the Pearl River Estuary (PRE) and adjacent seas. To provide valuable information to assist the governors, stakeholders, and decision-makers in tracking changes in environmental conditions, daily nowcasts and two-day forecasts from the ecological prediction system, namely the Coupled Great Bay Ecological Environmental Prediction System (CGEEPS), has been setup. These forecast systems have been built on the Coupled Ocean–Atmosphere–Wave–Sediment Transport modelling system. This comprises an atmospheric Weather Research Forecasting module and an oceanic Regional Ocean Modelling System module. Daily real-time nowcasts and 2-day forecasts of temperature, salinity, NO2 + NO3, chlorophyll, and pH are continuously available. Visualizations of the forecasts are available on a local website (http://www.gbaycarbontest.xyz:8008/). This paper describes the setup of the environmental forecasting system, evaluates model hindcast simulations from 2014 to 2018, and investigates downscaling and two-way coupling with the regional atmospheric model. The results shown that though CGEEPS lacks accuracy in predicting the absolute value for biological and biogeochemical environmental variables. It is quite informative to predict the spatio-temporal variability of ecological environmental changes associated with extreme weather events. Our study will benefit of developing real-time marine biogeochemical and ecosystem forecast tool for oceanic regions heavily impact by extreme weathers.


Introduction
The Pearl River Estuary (PRE) in South China is a large subtropical permanent open estuary fed by the Pearl River, which is the second largest river in China in terms of freshwater discharge (3.3 × 10 11 m 3 /yr) (Figure 1). Due to the substantial amount of anthropogenic nutrient inputs along with high discharge, the PRE is relatively productive with a high level of biodiversity. It provides crucial habitat and natural resources for multiple marine and freshwater fish species and the Chinese white dolphin, which is designated as a first-class national protected animal (Jefferson and Hung, 2004;Zhou et al., 2019). Oyster aquaculture has also developed rapidly in the PRE since the Ming and Qing dynasties, generating high revenue for the Guangdong government (Liu and Song, 2022). According to Cai and Li (2011), the PRE marine ecosystem is estimated to be valued at $31 billion per year. This represents an important component of the Gross Domestic Product of the Guangdong-Hong Kong-Marco Greater Bay Area.
Over the last few decades, rapid economic growth and progressive urbanization in the Pearl River Delta have led to considerable population growth in large cities in this region, including Hong Kong, Macao, Shenzhen, Zhuhai, and Guangzhou. This population increase poses a threat to the water quality of the PRE due to dumping of new pollutants, such as organophosphorus compounds, pharmaceuticals, and microplastics derived from personal care products (Stedman and Dahl, 2008;Zhou et al., 2016;Wang and Rainbow, 2020). These pollutants can cause severe damage to estuarine ecosystems and reduce biodiversity and fishery resources. This has caused substantial ecological and economic losses. Thus far, the loss value of the PRE is estimated to be approximately $5 billion per year. Therefore, protecting the estuarine environment from further water quality degradation is of considerable interest to stakeholders.
Real-time environmental forecasting systems based on threedimensional, coupled hydrodynamic-biogeochemistry models are important tools for providing information on spatio-temporal varied environmental conditions (Wang, 2001). These data can be used to assist local governments and private partners to establish specific strategies for water quality improvement. The increase in anthropogenic derived nutrients has fueled larger and more frequent hypoxia in recent years (Li D et al., 2020;Ma et al., 2009;Ye et al., 2017;Qian et al., 2018). Oyster aquaculture has been shown to mitigate eutrophication and hypoxia, whereas the location and size of the farming area require circulation and appropriate biogeochemical conditions (Yu and Gan, 2021).
The coupled hydrodynamic-biogeochemical model usually serves as the internal core of a real-time environmental forecasting system. Web-based visualization systems are another important component that aid in forecasting. Real-time environmental systems with visualization systems can provide an intuitive way for governments and private partners to understand data and facilitate rapid decision-making. These systems have been established in many estuarine areas and oceans adjacent to coastal areas worldwide. For example, the Chesapeake Bay Environmental Forecast System (CBEFS) runs every 6 h over 1-2 km as well as the Chesapeake Bay Regional Ocean Modelling System (ChesROMS) -Estuarine Carbon Biogeochemistry (ECB) Model (Feng et al., 2015). CBEFS can provide real-time nowcast and two-day forecasts of salinity, water temperature, pH, aragonite saturation state, alkalinity, dissolved oxygen, and hypoxic volume in graphics on the local website (www.vims.edu/hypoxia). Another example of a real-time forecasting system is the China high-resolution coastal ocean ecological environment numerical prediction system. This was constructed using the Regional Ocean Modelling System-Carbon Silicate and Nitrogen Ecosystem model (ROMS-CoSiNE) and Semiimplicit Cross-scale Hydroscience-Integrated System Model-CoSiNE (SCHIMS-CoSiNE) for the BoHuang Sea, Yellow Sea,South China Sea,and Xiamen Bay (http://202.108.199.24:8080/BJ_SZYB_Web/ default.html). The system can predict water temperature, current, multiple inorganic nutrients, chlorophyll, dissolved oxygen (DO), and pH predictions every 6 h. For the PRE, we constructed the Guangdong-Hong Kong-Marco Greater Bay Area Ecological Security Prediction System V1.0 using the Regional Ocean Modelling System (ROMS), Soil and Water Assessment Tool (SWAT), and Estuarine Carbon Biogeochemistry Model (ECB). The system was run every 24 h and primarily used for forecasting ocean currents, temperature, salinity, NO 2 + NO 3 , chlorophyll, and DO (http://zhujiangtest.xyz:8000/). This study describes an updated version of the real-time Ecological Environmental Forecast System for the Pearl River Estuary, the Coupled Great Bay Ecological Environment Prediction System (CGEEPS). We have introduced a general model configuration and evaluated the model accuracy for multiple hindcast variables and graphic visualization. CGEEPS provides daily nowcasts and two-day forecasts of environmental conditions throughout the PRE for temperature, salinity, NO 2 + NO 3 , chlorophyll, DO, and pH. Configuration of the first generation of CGEEPS used a coupled ocean circulationatmosphere model with a nitrogen-based ocean biogeochemistry component. The river inputs were generated from the land surface model SWAT climatology. We have added a range of updates that build on the first generation model (http://zhujiangtest.xyz:8000/). This includes a carbonate chemistry module that has been previously implemented in an existing nitrogen-based biogeochemistry model. River discharge was updated in near realtime from the China National Water and Rain Information Center (http://xxfb.mwr.cn/) and snatched by the Python request module.
The remainder of this paper is organized as follows. Section 2 presents the configuration of the modelling system core for the Forecast System and the model validation methodology. Section 3 presents the model validation and advantages of coupling it with the atmospheric model in a typhoon passage case. Section 4 summarizes the results and discusses the disadvantages of the current version of the forecasting system and potential future techniques that can be implemented.

Atmosphere-ocean circulation-carbon biogeochemistry model implementation
The physical model of CGEEPS is built on the Coupled Ocean-Atmosphere-Wave-Sediment Transport (COAWST) modelling system. This comprises an atmospheric Weather Research Forecasting (WRF) module and a ROMS (Warner et al., 2008;Warner et al., 2010). The WRF domain spans from 98°to 122°E and from 16°to 30°N with an average horizontal resolution of approximately 12 km. This encompasses the entire Pearl River Watershed and northern part of the South China Sea. The ROMS domain ranges from 105°to 121°E and from 16°to 26°N, covering the continental shelf of the North South China Sea (NSCS) and the PRE (Figure 2A). The horizontal grid spacing varies between 2.8 km and 40.7 km with an average resolution of approximately 4 km. The model has 30 vertical layers that follow the terrain with a higher resolution near the surface and bottom boundaries. The vertical scoordinate function is based on that of Shchepetkin and McWilliams (2009). The physical model is configured to use the recursive MPDATA 3-D advection scheme for tracers, four-order horizontal advection of tracers, third-order upwind advection of momentum, and the Mellor and Yamada (1982) turbulence closure scheme for vertical mixing. A viscosity/diffusivity coefficient of 10 -5 m 2 s -1 was used for biharmonic horizontal mixing in the momentum and tracer equations. Based on the CFL criterion, which is one of the global stability conditions, the internal mode time step of the numerical integration was 900 s, while the external mode time step was 60 s.
The circulation model was coupled with a nitrogen cycle model with consideration of carbonate chemistry and O 2 (Fennel et al., 2006(Fennel et al., , 2008(Fennel et al., , 2011. The model has 12 state variables, namely phytoplankton (Phy), chlorophyll (Chl), zooplankton (Zoo), nitrate (NO 3 − ), ammonium (NH 4 + ), O 2 , dissolved inorganic carbon (DIC), total alkalinity (TAlk), and two detritus pools. These include small and large detritus, each split into nitrogen (SDetN and LDetN) and carbon (SDetC and LDetC) ( Figure 2B). The mathematical description of the biogeochemical source and sink terms for the state variables is as follows.
Phy −n NH + 4 + l BM Zoo + l E Phy 2 k P + Phy 2 bZoo +r N SD SDetN + r N LD LDetN (10) ∂ DIC ∂ t = − q P C:N m max (T)f (E)L N Phy +r C SD SDetC +r C LD LDetC + r C RD RDOC + q Z C:N (l BM + l E Phy 2 k P + Phy 2 b)Zoo (16) The parameter definitions, values, and units used in this model for coastal applications are listed in Table S1. The inclusion of inorganic carbon and alkalinity as state variables is critical for the successful simulation of pH using CO2SYS (Lewis and Wallace, 1998), as implemented in MATLAB. Phosphate is not included in this version of the model but will be added in the future to improve the system.
The bathymetry of the model domain was obtained from ETOPO and smoothed to reduce truncation errors, with the minimum water depth set to 5 m. The minimum depth was set to 5 m to ensure global stability, i.e., h min +z max >0 (h min is the min undisturbed water depth and z max is the max free surface elevation at the coast due to gust wind-induced storm surges, Wang, 1996). The river discharge data were sourced from the three main branches of the Pearl River (Xijiang, Dongjiang, and Beijing Rivers) and from the Ministry of Water Resources of the People's Republic of China (http:// xxfb.mwr.cn). Freshwater inflows for the forecast period of each daily simulation were held constant for the two-day forecast. This was based on the inflow during the nowcast period. The data were summed and assumed to be distributed in the upper part of Lingding Bay, the lower Modaomen sub-estuary, and Huangmaohai Bay. River nutrient concentration was retrieved from water samples collected monthly from the eight river outlets of the Pearl River Delta (Lu et al., 2009) This included NO 3 , NH 4 , and organic nitrogen. All riverine organic matter was treated as SDetN in the simulations. The phytoplankton, zooplankton, and LDetN values were set to zero. DIC and alkalinity riverine inputs were sourced from the Xiamen University on request. Given that the nutrient and carbon-associated measurements are highly limited, climatological values were used at the open boundary, and the physical module was forced by the daily average temperature, salinity, current velocity, and de-tidal water level from a high-resolution regional South China Sea model (Peng et al., 2015). The biochemical variables used for the boundary conditions of the model were derived from the historical simulation of the CMIP6 model CESM-WACCM. The main inputs and open boundary conditions are presented in Table 1. We spun up the model for over six years before using the available historical cruise data (2014-2018) to validate the model.

Forecast system configuration
CGEEPS is maintained manually and runs at a specific time of the day to automate the workflow of the prediction system. The forecast system retrieves the necessary information from NOAA ftp to download the seasonal climate forecast from the NCEP coupled forecast system model Version 2 (CFSv2) using the "wget" command. Python was used to interpolate the downloaded data to a fine WRF atmospheric model grid. The Python requests module was used to snatch the river discharge from the China National Water and Rain Information Center. NetCDF operators combined with MATLAB and Python were used to modify the NetCDF input and output files efficiently. The forecast system maintainer replaces dates in a generic ROMS text input file for the simulation to commence on the correct day. The system was run on a workstation with 66 CPUs and simulated for three days, including a one-day nowcast and a twoday forecast. Each successive nowcast resulted from the end of the nowcast for the previous day. MATLAB was used to further postprocess the model NetCDF output to generate portable network graphics (png) files for online visualization.
An ecological security prediction system requires nowcast and forecast boundary conditions. The WRF model was used to produce real-time forecasting meteorological conditions for CGEEPS ( Figure 3). The WRF provides evaporation and precipitation, sea surface heat flux, sea surface wind, and air pressure data to the ocean model. The ROMS subsequently feeds the sea surface temperature back to the meteorological model. The exchange frequency was once every 10 min. In this study, the impact on the study area of Hurricane Hato in August 2018 was simulated with and without the WRF.

Observational data
We used a series of cruise observations from 2014 to 2018 to validate the model ( Figure 4; Table 2). The observations are part of an ongoing program from a ship time-sharing project sponsored by the National Natural Science Foundation of China (NSCF) and Xiamen University. Measurements from Conductivity, Temperature, and Depth (CTD) sensors documented the vertical structure and monthly temperature, salinity, and DO. Chlorophyll, NO 3 , NH 4 , and pH levels were measured in the laboratory.

Model skill metrics
Quantitative model-data comparisons using multiple skill metrics are critical because they highlight the advantages and potential limitations of a particular model, which must be carefully considered before using the model as a tool for scientific study or decision-making Stow et al., 2009). In this analysis, multiple skill metrics were examined. This included the correlation coefficient the that measures the tendency of the model and observed values to covary, and was calculated as O is the mean of the observations, M is the mean of the model estimates, and n is the total number of observations available for comparison with model estimates. The bias, unbiased root-meansquared difference (unbiased RMSD), and total root-mean-squared difference (RMSD) were calculated as: These three skill assessment statistics are particularly useful as bias reports of the size of the model-observation discrepancies. Bias values near zero indicate a close match. However, this can be misleading because negative and positive discrepancies can cancel Flow chart of the forecast system automated workflow. each other. The unbiased RMSD nullifies the effect of the mean and is a pure measure of how model variability differs from observational variability. Total RMSD provides an overall skill metric, given that it includes components for assessing mean (bias) and variability (unbiased RMSD).

Target and Taylor diagrams for model validation
Multiple model skill metrics are compactly visualized on Taylor (Taylor, 2001) and target diagrams (Hofmann et al., 2008;Friedrichs et al., 2009;Jolliff et al., 2009). For the Taylor diagram, in addition to these skill metrics, the model standard deviation, s m has also been calculated.
For both types of diagrams, skill statistics are typically normalized by the observational standard deviation (s o ). This allows plotting of multiple different datasets on the same diagram.
The Taylor diagram is plotted using a polar coordinate system and summarizes the skill metrics, including r, s, and unbiased RMSD. In contrast, the target diagram was plotted using a Cartesian coordinate system and summarized the total RMSD, unbiased RMSD, and bias. On the normalized Taylor diagram, the reference point at (1, 0) represents a perfect skill score. Meanwhile, on the target diagram, the center of the target (0, 0) represents a perfect skill score.

Hindcast simulations for evaluating model accuracy
The hindcast simulation of sea surface salinity was in strong agreement with the observations for all eight cruises ( Figure 5). The correlation coefficient was 0.89. The bias and RMSD were 0.83 and 2.77 PSU, respectively. The simulated sea surface salinity showed Comparison between observed (dots) and simulated (background color) surface salinity (< 2 m) for datasets (1) to (8)   . The simulated plume morphology is consistent with previous observational facts reported in literature that the Pearl River buoyant plumes formed during summer can be classified into four types, namely offshore bulge spreading, west alongshore spreading, east offshore spreading, and symmetrical alongshore spreading. These were predominantly determined by the combined effect of river discharge and wind (Chen et al., 2017;Luo et al., 2012;Zu and Gan, 2015;Ou et al., 2009). River discharge has regulated the plume size, while wind has played an important role in changing the plume shape. The east and southeast winds drive the buoyant plume westward, resulting in the plume spreading westward alongshore. Given that the river plume resides at the surface, it is susceptible to the surface Ekman effect. The southwest wind upwelled and detached the eastward plume off the coast, which then spread offshore. If the wind is purely in a southerly direction, the plume is confined nearshore and moves both eastward and westward, thus forming a symmetrical structure. Next, we compared the CGEEPS hindcast temperature and salinity along the vertical transect for six cruises where profiled CTD data were available ( Figure 6). The simulated temperature and salinity were in strong agreement with the observations. For temperature, the CGEEPS hindcast replicated the observational fact that there was a strong thermocline in September 2015 and 2016, June 2017 and 2018, and August 2019, but not in October 2017. The correlation coefficient was 0.81. The bias and RMSD was 0.36 and 1.74°C, respectively. The formation of a strong thermocline was due to upwelling of cold water under southwesterly winds, which generally occurs from June to August.
In addition to the thermocline from upwelling, the low salinity riverine water extended further offshore with the southwesterly wind in both the CGEEPS hindcast simulation and observations. The 28 psu isohaline extended over 22°N. However, October 2018 was highly favorable for the occurrence of downwelling conditions. The low salinity water was increasingly confined adjacent to the river mouth and the isohalines were vertically distributed. The correlation coefficient was 0.94 and the bias and RMSD were 0.19 and 1.78 psu, respectively.
The high level of skill in terms of temperature and salinity fields suggested that our modelling system was relatively robust in replicating the dynamics of the PRE-ocean system, which is a river-dominated margin (Rabouille et al., 2008).
We subsequently examined how the nitrogen type nutrients, which are a combination of NO 3 + NO 2 + NH 4 , were reproduced by CGEEPS. Given that rivers and upwelling are the two main sources of nutrients, we initially observed that the concentrations of N-type nutrients can reach as high as 150 mM inside the estuary and decrease to less than 1 mM as the plume water extends offshore. The nitrate concentration was relatively high (5-10 mM) at 20 m from the bottom owing to upwelling (Figure 7). The correlation coefficient between simulated and observed NO 3 was 0.91; however, the bias seemed to be slightly larger (Bias = 20.79 mM, RMSD = 63.1 mM). We examined the output data and found that NO 3 -overestimation predominantly occurred at three stations in the upper reaches of the estuary closer to the river mouth boundary of the model grid. This overestimation was likely because we had set up at the river mouth in the upper part of Lingding Bay rather than at the actual river mouth. If we excluded these three points from the calculation, the bias and RMSD were almost halved (Table 3).
The nutrient distribution was in line with the salinity distribution pattern in the river plume region, which were the two variables with the best performance in the simulation system.
The full model observation skill assessment is presented in Table 3. The comparison was conducted by gathering all the observational points and extracting the model output at the observational stations. This implies that the skill number incorporated both spatial and temporal variability. In addition to the temperature, salinity, and N-type nutrients, we also examined the model performance for surface chlorophyll concentration and bottom DO (Table 3). The comparison results showed that the correlation between the observed and simulated surface chlorophyll content was 0.29. The surface chlorophyll concentration was overestimated by approximately 3 mg-Chl/L, particularly in the plume area. This was likely due to the deficiency of the current ecosystem module in CGEEPS. To date, we have not incorporated the P-and light-limitation effects due to resuspended sediment. We expected that further incorporation of these processes would slow the growth of phytoplankton and decrease the surface chlorophyll concentration. Comparison between observed (dots) and simulated (background color) temperature (A-F) and salinity (G-L) along vertical section for datasets (1) to (6).
Although the surface chlorophyll concentration had been overestimated, the subsurface high chlorophyll maximum (DCM) structure was reproduced by the model (not shown here). The DCM layer lies immediately above the upwelled high NO 3 concentration, where nutrients and light were determined to be sufficient.
We further estimated the model's overall ability to reproduce bottom DO and surface pH. The model's skill in producing surface chlorophyll concentration and bottom DO was relatively poor. The correlation coefficient was positive at 0.41 for surface pH, but negative at −0.41 for bottom DO. The surface pH was underestimated by approximately −0.11 and the bottom DO was overestimated by approximately 1.2 mg/L. The underestimation of pH was likely associated with the surface temperature being underestimated by the model. This resulted in more dissolved CO 2 being held in the water and more hydrogen ions being released, resulting in lower pH values. The bottom DO was overestimated by the model, which was unexpected because the surface DO chlorophyll overestimated and more oxygenconsuming materials should be available. This overestimation may be because the current version of the biogeochemistry model did not include sediment oxygen consumption (SOC). Currently, the biogeochemistry model assumes that organic matter is immediately consumed upon reaching the bottom, whereas measurements show that organic matter consumption can be taken up over months. Instead of remaining in the same location, it is transported as fluid mud and undergoes numerous cycles of resuspension and deposition. The movement and long-term retention of organic matter in the sediment layer consumes additional oxygen when decomposed by bacteria. The current version of the carbon biogeochemistry model does not include the dissolved form of organic matter, which is another significant pool reserving organic matter and consuming DO. The current model parameter selections are predominantly empirical without optimization, which introduces considerable uncertainty to the simulated biogeochemical variables. Table 3 provides the overall skill of the model by compiling multiple cruise data. However, it was not known how the model reflects spatial and temporal variability. Next, we performed a cruiseby-cruise comparison, which provided more specific information on the model's skill in replicating the spatial variability of multiple variables. Both Taylor and target diagrams were introduced to graphically and quantitatively visualize multiple model skill metrics (Figure 8). For the predicted environmental variables at the surface and on the bottom, the model performances for salinity and Comparison between observed and simulated nitrogen type nutrients (NO 3 + NO 2 + NH 4 ) for the surface (A-D) and along vertical transect (E-H) for datasets (3) to (6).
temperature were uniform for all cruises. For salinity, the bias and RMSD were lower than 2 and 3.6 PSU, respectively, but the correlation coefficients were higher than 0.9. However, for temperature, both the bias and RMSD were low, and the correlation coefficients were also low ( Figures 8A, B). In contrast to salinity and temperature, the performance of biogeochemical variables, especially surface chlorophyll, was cruise dependent. The model output showed a small bias and RMSD in October 2017, but a large bias and RMSD in June 2017. Similarly, the correlation coefficient was −0.08 in October 2017 but reached 0.61 in August 2018. Although the overall correlation coefficients for surface pH and chlorophyll concentration were 0.41 and 0.29, respectively, the correlation coefficient for a single cruise can reach as high as 0.7. These results suggest that surface chlorophyll and pH are more challenging to predict than temperature and salinity, because ecological and biogeochemical processes generally have larger uncertainties than physical processes. Sectional data are only available for salinity and temperature. They showed a small bias and RMSD for all cruises. The correlation coefficients are over 0.7 for most cruises, except in September 2016, when the correlation coefficient for temperature was approximately 0.4.
Currently, data from multiple cruises are not collected at the same locations. Therefore, the ability of the model to reproduce temporal variability has not yet been assessed. Raw cruise data are processed using multiple laboratories. This implied that we could not guarantee that the same criteria and standards would be followed during data collection and preliminary processing. Therefore, the data quality itself may degrade the overall performance. We expect a future observational network and data-sharing platform to be set up in this region, which would benefit the systematic evaluation of the simulation performance in the future.
3.2 Differences in physical and biogeochemical fields by coupling with regional atmospheric model Most coupled ocean-biogeochemistry model use atmospheric reanalysis products as forcing. One advantage of CGEEPS is the coupling of the regional atmospheric model with the regional ocean model. Compared to reanalysis products, the regional model usually has a high spatio-temporal resolution and the feedback of SST to the atmosphere. The coupled model has unique advantages in terms of simulating hurricane tracks and intensity. Therefore, we examined how the coupling technique has impacted ocean dynamics and biogeochemistry by using Hurricane Hato (2017)  landfall, Hato weakened rapidly to 30 ms −1 over several hours and dissipated by 0800 UTC on August 24.
Two simulations were conducted from August 07 to September 05, 2017, the month when Hurricane Hato (2017) passed through. One was undertaken using two-way coupling with high-resolution regional configured WRF and the other by one-way forcing with coarse-resolution global ERA5. The simulated tracks of typhoons across the NSCS during August 22-24, 2017 are shown in Figure 9. The differences in both physical and biogeochemical fields are shown in Figures 10, 11.
We initially found that the typhoon tracks under the two runs were comparable. This indicates that the track of the typhoon was predominantly controlled by the background of large-scale atmospheric circulation (Figure 9).
For both runs, SST decreased and SSS increased rapidly from August 22 to August 23, 2017 ( Figures 10B, C). This was caused by strong vertical mixing induced by the high wind of the hurricane. The SST increased after August 23 and decreased again from August 26. This was caused by upwelling after the eye passed offshore. The magnitude of the temperature drop from August 26 to August 28, 2017, in the two-way coupling WRF run was smaller than that in the one-way ERA5 run. This was likely because SST feedback to the atmosphere during the coupling weakened the strength of the tropical cyclone. The SST recovered from 28 August in both runs. However, the coupled WRF run diverged from the ERA5 one-way run from 31 August by maintaining a high SST at approximately 29°C. At the same time, the tropical cyclone brought heavy precipitation. This resulted in lower SSS in the twoway coupled WRF run than in the ERA5 one-way run.
Due to SST feedback to the atmosphere, we found that the difference between the two runs for surface air temperature (0.57°C ) was larger than SST (0.34°C) ( Figure 10A). The sea surface salinity was maintained at approximately 32°C until August 27, 2017, and subsequently decreased. This was because river runoff increased after hurricane-induced heavy rainfall, which introduced a substantial amount of freshwater to the region being studied. The simulated results are consistent with the previous studies on ocean SST cooling and enhanced vertical mixing during a passage of a typhoon (Shen et al., 2021).
Surface chlorophyll, bottom DO, and surface pH varied with SSS and SST (Figures 10D-F). The temporal variability of chlorophyll, DO, and pH are strongly related to each other, that is, higher chlorophyll concentrations result in lower DO and higher pH. High surface chlorophyll concentrations indicate high primary production, which indicates that more oxygen-consuming materials are generated and DO decreases. A higher primary production also results in a higher pH because the equilibrium of carbonate ionization moves left after CO 2 is taken up by photosynthesis, thereby decreasing the concentration of hydrogen ions and increasing the pH value.
Owing to the strong vertical mixing induced by the hurricane, the surface chlorophyll concentration decreased, bottom DO increased, and surface pH decreased (Figures 10D-F). During the upwelling period when SST decreased, surface chlorophyll increased, bottom DO decreased, and surface pH also increased, which was caused by high primary production stimulated by subsurface high nutrients. When the SSS decreased, the surface chlorophyll increased, the bottom DO decreased, and the increase in surface pH became larger.
A comparison between the two-time series and the two runs showed that the surface chlorophyll concentrations yielded larger differences than the SST and SSS differences. During Hato's (2017) passage, although the chlorophyll concentration decreased in both runs, the average surface chlorophyll concentration in the ERA5 run was approximately 2.5 mg-Chl/L higher. The concentration differences were approximately 6.0 mg-Chl/L two weeks later (September 02, 2017). This suggests that the physical and biological processes are highly nonlinear. A minor change in the physical field may be enlarged several times in biogeochemical fields.
For a better view of the differences, we plotted the horizontal distribution of these variables for August 23 and September 02, The observed and simulated track of Hurricane Hato (2017). The six-hour best-track TC data were obtained from the Shanghai Typhoon Institute.
2017. During hurricane Hato's (2017) passage and 11 d after it had passed, was when increased runoff substantially impacted biogeochemical processes.
On August 23, 2017, typhoons Hato induced strong vertical mixing and onshore currents, resulting in low-temperature water mixing up to the surface and confinement to the narrow nearshore band (Figure 11). Low salinity and riverine waters were predominantly confined to Lingding Bay. As a result, there was relatively higher chlorophyll and lower DO within Lingding Bay and adjacent to the nearshore western bank. As river discharge increased and was delivered to the coastal ocean, the plume water extended offshore. As a result, high surface chlorophyll, low DO, and high pH water also extended further offshore. The low salinity, high chlorophyll, and bottom DO water in the one-way ERA5 run occupied a larger area than that in the coupling WRF run. This was likely because the regionally configured WRF had a higher resolution than the one-way ERA5, which induced high sub-grid vertical mixing.

Forecasts for eutrophication, hypoxia, and ocean acidification
Although our simulation system may not have an exact accuracy in replicating the absolute value of all environmental variables at present, it is relatively informative in predicting variability under extreme weather conditions, such as floods and droughts.
The Pearl River Basin experienced a rainstorm in June 2022, which rapidly increased river discharge to 6×10 4 m 3 s −1 on June 15 ( Figure 12A). CGEEPS successfully captured the surface salinity decrease associated with a sudden increase in runoff ( Figure 12B). Compared to June 01, 2022, a low-salinity bulge formed on June 15 beside Lingding Bay ( Figure 13B). The sudden increase in runoff also delivered more nutrients from land to the estuarine shelf, leading to phytoplankton blooms and a decrease in bottom dissolved oxygen ( Figures 12E, F). At the surface, the low-salinity bulge was accompanied by a high NO 3 concentration ( Figure 13D). However, the surface chlorophyll water was not synchronized with surface salinity and NO 3 , which was relatively low within Lingding Bay and nearshore, but relatively high offshore ( Figure 13H).The surface chlorophyll concentration depends on the relative importance of flushing and phytoplankton growth. High runoff reduces the residence time of phytoplankton, flushing them nearshore. Although the nutrient concentration was relatively high, there was insufficient time for phytoplankton to take up nutrients. Therefore, surface chlorophyll concentration was maintained at a low level. In contrast, phytoplankton had sufficient time to uptake NO 3 and bloom in the offshore region. We observed that high surface chlorophyll water was accompanied by higher pH levels. The pH was regulated by primary production. River runoff also brought a mass of DIC, resulting in a significant decrease in pH adjacent to the river mouth ( Figure 13H). We also observed a sudden decrease in pH on June 16 from the time series after the peak discharge ( Figure 12D). These results suggest that the PRE and adjacent seas became more eutrophic, hypoxic, and acidic due to flooding. CGEEPS can capture such temporal and spatially varying environmental conditions, which can facilitate stakeholders making daily decisions regarding their usage of the bay's resources.
To make it easier for the governments and private users to access and understand relevant information, the information provided by the CGEEPS system should be presented in an easily accessible format. In addition to the near-real-time simulation system, CGEEPS also includes a real-time online data visualization system in which the visualized graph is published in a webpage format. Graphics generated from CGEEPS can nowcast environmental variables daily and forecast them over two days (Figures 14, 15). The focus has been on the horizontal distribution of variables in the first generation of the forecast system (http:// zhujiangtest.xyz:8000/). To better inform governance and stakeholders, we have recently extended it to produce a vertical profiles and maps to display the subsurface structure, which is informative for visualizing eutrophic, hypoxic, and acidic water.

Discussion
Daily nowcasts and two-day forecasts from the ecological prediction system, the CGEEPS, were set up in the PRE and adjacent seas. The forecast systems were built on the COAWST modelling system. This comprises an atmospheric WRF module, an oceanic ROMS module, and a carbon-based biogeochemistry module (FENNEL). The forecast system can now provide realtime nowcast and two-day forecasts for temperature, salinity, NH 4 + NO 2 + NO 3 , chlorophyll, DO, and pH. Hindcast comparison with available historical data from summer 2015-2018 has shown that the prediction system can replicate the temperature, salinity, and Ntype nutrients for surface and vertical transects. However, the model performance for chlorophyll, DO, and pH was relatively poor. We also compared how the predicted biogeochemical fields differed from those predicted using coarse-resolution atmospheric reanalysis products. Using Hurricane Hato (2017) as an example, we found that stronger sub-grid vertical mixing was introduced after coupling with WRF. As a result, the surface salinity was lower than that when using coarse-resolution one-way ERA5 forcing. Due to the non-linear interaction between the physical and biogeochemical fields, the differences in surface chlorophyll became much greater than the salinity. Results suggest that the ecological environment prediction system was highly sensitive to physical forcing. The improved physical model would be beneficial to the system accuracy. We also determined how the physical and biological variables were predicted during the June 2022 rainstorm event. Results have shown that CGEEPS successfully captured data showing that the PRE and adjacent seas became fresher, more eutrophic, hypoxic, and acidic due to the Pearl River floods associated with the rainstorm. Results suggest that CGEEPS is capable of capturing the spatio-temporal variability of ecological environmental changes associated with extreme weather events. At present, our ecological prediction system lacks accuracy in predicting some biological and biogeochemical environmental variables. We expect future improvements to increase system predictability, including: 1) more precise description of the water column and sediment processes, and 2) optimizing the parameter scheme of biogeochemical models. In addition to the modelling system itself, we found that building an observational network and data-sharing platform for this region was important. Compared to the Texas-Louisiana shelf, which is heavily impacted by the Mississippi River and Chesapeake Bay in the US, current observational data used to validate the system are relatively sparse and not spatially uniform. This increases the difficulty in evaluating the performance of the prediction system after new processes are introduced. An observational network has served as the basis for better parameterization of the ecosystem model. Another disadvantage of CGEEPS is that the system is currently operated and maintained manually, and we expect to use the Cron software utility to completely automate the CGEEPS workflow in the next phase of development. Despite these disadvantages, CGEEPS is the first ecological environmental prediction system for the China Great Bay Area. The Hong Kong University of Science and Technology model data platform (http://ocean.ust.hk:8080/SiteMapApi/new/index.jsp) provides data visualization for daily averaged ocean circulation, ecosystem, and hypoxia, but no forecasts for these fields. Another example is the real-time Regional Forecast System of the SCS Marine Environment (https://epanf.scsio.ac.cn/), which is used to predict ocean circulation, waves, and several other important oceanatmospheric dynamic variables rather than ecological environmental conditions (Zhu et al., 2020). Real-time environmental forecasting systems based on coupled oceanatmosphere have also been set up worldwide. The Coupled Northwest Atlantic Prediction System (CNAPS), which is run by the North Carolina State University, predicts marine weather, ocean waves, and ocean circulation daily over an extensive area of the coastal northwest Atlantic Ocean. Another example is the South Eastern Levantine Israeli Prediction System (SLEIPS), which was setup through the Princeton Ocean Model (POM) and Variational Initial and Forcing Platform (VIFOP) for the coastal zone of the south eastern corner of the Levantine Basin (https:// isramar.ocean.org.il/isramar2009/selips/), which provides temperature, salinity, and sea current prediction every 24 h. Compared with these prediction systems, CGEEPS is capable of predicting ecological environmental variables based on the coupled ocean-atmosphere-biogeochemistry technique. However, many Real-time graphics displaying the surface (A) and vertical transect (B) of salinity (Screen shot from July 07, 2022). challenges on prediction accuracy were present because biogeochemical models are highly unconstrained and observation streams are much sparser than physical streams (Fennel et al., 2019).
Therefore, we implemented a coupled ocean-atmosphere technique for ecosystem and biogeochemistry prediction. However, for river-dominated ocean margins, such as the Pearl River estuarine-coastal system, the land delivers substantial amounts of nutrients, organic matter, and suspended sediments. These materials regulate the biogeochemical cycle of estuarinecoastal ocean systems; therefore, coupling with process-based land ecosystem models is required in the future. Such a land-estuaryocean biogeochemistry model has been developed for the Chesapeake Bay system (Feng et al., 2015). However, only land delivered materials in one-year was cycling used, and river discharge was still the USGS measurements (Bever et al., 2021). Model coupling techniques, such as CPL7, are required to develop an operational system (Sun et al., 2020).
Coastal marine ecosystems are subject to multiple stressors, such as climate change and human land activities. Implementation of an ecological environmental prediction system is important for protective and adaptive measures in ocean ecosystems. The prediction tool that we have constructed is based on realistic representations of ocean circulation coupled with biogeochemical and ecological models, which can forecast short-term (days to weeks) to seasonal (months) time intervals. In addition to this type of model, statistics based artificial intelligence models have emerged and have been applied to coastal ocean biogeochemistry studies in recent years (Chen et al., 2019;Li X et al., 2020;Ou et al., 2022;Yu et al., 2022). Compared with the coupled hydrodynamicbiogeochemistry model, artificial intelligent models are more efficient and require less computational resources. However, the application of such models to ocean biogeochemistry is still in its early stages. Therefore, the predictability and results explainable for such a model require further investigation. Real-time graphics displaying the surface (A) and vertical transect (B) of nitrate (Screen shot from July 07, 2022).

Concluding remarks
Warming, acidification, deoxygenation, and eutrophication of estuarine-coastal oceans have manifested in recent decades due to climate change and human land activities. Ecological environmental prediction systems are important tools that help decision-makers and the public implement protective and adaptive measures for ocean ecosystems. Based on the coupled ocean-atmospherebiogeochemistry model, we constructed a real-time ecological environmental forecast system for the Pearl River Estuary (PRE), namely the Coupled Great Bay Ecological Environment Prediction System (CGEEPS). Compared to previous ecological prediction systems built for this region, CGEEPS can not only provide realtime nowcasts and two-day forecasts for chlorophyll but also DO and pH. Although limited in its absolute accuracy, the system is still informative for eutrophication, deoxygenation, and acidification, especially under extreme weather events. We expect the expansion of biogeochemical and ecological observational systems built for this region to help develop and apply such prediction tools in the future. Ultimately, the predicted data products can benefit the economic, environmental, and public safety needs of stakeholders and governments.

Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below:; http://www.gbaycarbontest.xyz:8008/.

Author contributions
LL, ZM: These authors contributed equally to this work and share first authorship. They are responsible for the writing of this paper. WM, JH, YLi, YLiu, YH, YZhu: These authors contributed equally to this work. They are responsible for the setup, calibration, operation and maintenance of the Coupled Great Bay Ecological Environmental Prediction System. YF: She share senior authorship. She is in charge of the whole project. All authors contributed to the article and approved the submitted version.