## 2 Materials and Methods
### Hydrodynamic Model  
The transport model was based on the currents calculated by a hydrodynamic model. Given that the coverage area of brown tides is restricted to the shallow intertidal and sub-tidal zones, a 2 DH model mesh was adequate to simulate the hydrodynamics of the study area, considering the balance of computational speed and model precision in Qinhuangdao [42; 43]. The model solved the two-dimensional incompressible Reynolds-averaged Navier-Stokes equations with assumptions of Boussinesq and hydrostatic pressure. The computational mesh used in this research was unstructured triangular mesh with the spatial discretization adopting the cell-centered finite volume method. The continuity and momentum equations were converted into nonlinear shallow water equations by integrating them over depth, and the 2 DH equations are shown as follows:  
\[\frac{\partial h}{\partial t}+\frac{\partial hu}{\partial x}+\frac{\partial hv }{\partial y}=hS \tag{1}\]  
\[\frac{\partial hu}{\partial t}+\frac{\partial lu^{2}}{\partial x}+\frac{ \partial lvu}{\partial y}=fvh-gh\frac{\partial l}{\partial x}-\frac{h}{\rho_{ 0}}\frac{\partial p_{a}}{\partial x}-\frac{gh^{2}}{2\rho_{0}}\frac{\partial \rho}{\partial x}+\frac{\tau_{xx}}{\rho_{0}}-\frac{\tau_{bx}}{\rho_{0}}-F_{u} +hu_{s}S \tag{2}\]  
Figure 2: The nested model domain and mesh. (**a**) The regional model domain and mesh of the Bohai Sea; (**b**) The local model domain and mesh of Qinhuangdao. \[\frac{\partial hv}{\partial t}+\frac{\partial huv}{\partial x}+\frac{\partial hv^{2} }{\partial y}=-fuh-gh\frac{\partial\eta}{\partial y}-\frac{h}{\rho_{0}}\frac{ \partial p_{a}}{\partial y}-\frac{gh^{2}}{2\rho_{0}}\frac{\partial\rho}{ \partial y}+\frac{\tau_{sy}}{\rho_{0}}-\frac{\tau_{by}}{\rho_{0}}-F_{v}+hv_{s}S \tag{3}\]  
where \(t\) is the time; \(x\) and \(y\) are the Cartesian coordinates; \(u\) and \(v\) are the depth-averaged velocity components in the \(x\) and \(y\) directions; \(h\) is the sum of the still water depth, \(d\), and the time-varying water level, \(\eta\); \(S\) is the discharge of the point sources; \(f\) is the Coriolis parameter related to the angular velocity of Earth's rotation and the geographic latitude; \(g\) is the gravitational acceleration; \(\rho\) is the water density; \(\rho_{0}\) is the reference water density; \(P_{a}\) is the atmospheric pressure; \(\tau_{sx}\) and \(\tau_{sy}\) are the surface stresses induced by wind; \(\tau_{bx}\) and \(\tau_{by}\) are the bottom stresses; \(F_{u}\) and \(F_{v}\) are the horizontal stress terms due to eddy viscosity; and \(u_{s}\) as well as \(v_{s}\) are the discharge velocities of the point sources. The model adopted wet-dry dynamic boundary processing technology, and the critical water depths of dry and wet points were 0.005 m and 0.05 m, respectively. The computational time step was in the range of 0.1-30 s and was adjusted to meet the Courant-Friedrich-Lewy (CFL) condition with a limited value of 0.8 in the model calculation process, where \(CFL=\left(\sqrt{gh}+\left|u\right|\right)\frac{\Delta t}{\Delta x}+\left(\sqrt{ gh}+\left|v\right|\right)\frac{\Delta t}{\Delta y}\). The horizontal eddy viscosity coefficient used the Samagorinsky model to calculate the grid scale, and the Samagorinsky coefficient was 0.28. The initial velocity and water level were 0 m/s and 0 m, respectively. The Manning number was 74 m\({}^{1/3}\)/s, by which the model bed resistance was determined. Based on the Bohai Sea model, the boundary condition at the tidal inlet and open sea boundaries was adopted for the medium tidal range and Flather condition, respectively [42]. The boundary condition at the river outlet was set as the specified discharge value of 1 m\({}^{3}\)/s, and the initial condition was set as 0, including values of surface elevation, u-velocity, and v-velocity. The time step interval could self-regulate itself, ranging from 0.01 to 30 s to meet the CFL restriction as velocity changed with time at all computational nodes and time steps.