## 2 Materials and Methods
### Model  
The ocean model utilized in this study was the FVCOM model originally developed by the University of Massachusetts Dartmouth (UMASSD) and improved by the efforts of the Woods Hole Oceanographic Institution (WHOI) [27; 28; 29]. This model adopts an unstructured triangular mesh to solve the primitive equation through the finite-volume discrete method, which not only retains the high efficiency of finite difference calculation but also incorporates the flexibility of the finite element method to fit the shoreline ([[http://www.fvcom.smast.umassd.edu/](http://www.fvcom.smast.umassd.edu/)]([http://www.fvcom.smast.umassd.edu/](http://www.fvcom.smast.umassd.edu/)), accessed on 8 March 2022). In addition, time integration of either a mode-split solver or a semi-implicit solver is provided to be chosen. The vertical eddy viscosity and thermal diffusion coefficients are calculated using the General Turbulence Model (GOTM) in the vertical and the Smagorinsky turbulent parameterization in the horizontal [30]. The model grid established in this study, the same as that in [31; 32], covered the Pearl River Delta network, the Pearl River Estuary, and the northern South China Sea within a 100 m isobath range, with a minimum grid size of 10 m within the PRE and 2 km resolution in the radar observing area over the continental shelf (Figure 2). The topography data are extracted by Google Earth. The bathymetry data used in the model is acquired from Nautical Chart and bed sweeping data provided by CCCC (China Communications Construction Company) first harbor engineering company limited. In the vertical dimension, a terrain tracking coordinate system combining \(\sigma\) and S coordinates was divided into 20 layers, with resolutions up to 0.25 m. The model bottom friction stress is controlled by \(\left(t_{bx},t_{by}\right)=C_{d}\sqrt{u^{2}+v^{2}}(u,v)\), where \(t_{bx},t_{by}\) are the x and y components of bottom stresses. The drag coefficient \(C_{d}\) is determined by matching a logarithmic bottom layer to the model at a height \(z_{ab}\) above the bottom, i.e.,  
\[C_{d}=max\left(k^{2}/\ln\left(\frac{z_{ab}}{z_{o}}\right)^{2},0.0025\right) \tag{1}\]  
where \(k=0.4\) is the von Karman constant and \(z_{o}\) is the bottom roughness parameter. The length scale of \(z_{o}\) set in our model is 0.001. The river discharges of seven rivers, namely, Xijiang, Tanjiang, Beijing, Liuxi, Zengjiang, Dongjiang, and Hanjiang, were considered in the model. The climatic river discharge data in the model were based on the daily discharges at the Gaoyao, Shijiao, and Boluo hydrological stations during the 2003-2007 dry season (November to April) provided by the Guangdong Hydrological Bureau. The TPX08 tidal model, developed by Oregon State University, was utilized as the tidal-driven force at the open boundary [33]. The elevations along the open boundary were calculated by the Tide Model Driver (TMD) using the TPX08 tidal model. Eight major tidal constituents (M2, N2, S2, K2, K1, O1, P1, and Q1) were included in the tidal drive, which was proven to be reasonable for simulating tidal currents around the PRE [2; 8]. Hourly wind fields, as the wind drive, were obtained from Climate Forecast System Version 2 (CFSV2) of the National Centers for Environmental Prediction (NCEP) [34]. The initial fields of temperature and salinity were generated by the climatology from Simple Ocean Data Assimilation (SODA) ([[https://climatedataguide.ucar.edu/climate-data/soda-simple-ocean-data-assimilation/](https://climatedataguide.ucar.edu/climate-data/soda-simple-ocean-data-assimilation/)]([https://climatedataguide.ucar.edu/climate-data/soda-simple-ocean-data-assimilation/](https://climatedataguide.ucar.edu/climate-data/soda-simple-ocean-data-assimilation/)), accessed on 1 January 2022). The model was run for 2 months to stabilize the river discharges and model states. In the observation-model comparisons, the control experiment was performed by applying tidal drive, river discharge, and wind forces.