## 2 Materials and Methods
### Numerical Simulation Method  
#### 2.3.1 Environment Settings  
Before numerical simulation, the radius and initial velocity of the hydrothermal vent in the simulated environment need to be determined. According to the theoretical model of the hydrothermal plume (Equation (10)), the relationship between the buoyancy frequency and the buoyancy flux of the vent can be obtained [9]. \[Z_{neutral}=2\left(\frac{B_{exit}}{N^{3}}\right)^{1/4}, \tag{10}\]  
where \(Z_{neutral}\) is the neutrally buoyant plume height, \(B_{exit}\) is the buoyancy flux of the vent orifice, and \(N\) is the buoyancy frequency. Based on buoyancy frequency (Equation (11)), buoyancy flux (Equation (12)), temperature gradient value, density gradient value and previous studies, the initial velocity of the vent orifice is calculated. After sorting out the equation, the relationship between velocity and radius can be obtained (Equation (13)) [6]. \[N=\sqrt{\frac{g}{\rho(T_{bottom})}\cdot\frac{\partial\rho(z)}{\partial z}}= \sqrt{\frac{g}{\rho(T_{bottom})}\cdot\frac{\partial\rho}{\partial T}\cdot \frac{\partial T}{\partial z}}, \tag{11}\]  
Figure 3: CTD 05-2 station information. (**a**) Turbidity and potential temperature of CTD 05-2; (**b**) The density and specific heat capacity of seawater used in the numerical model change with temperature. Figure 2: Water anomaly detection in the Longqi hydrothermal field. (**a**) The turbidity outlier; (**b**) the potential temperature. \[B_{exit}=g\frac{\rho(T_{bottom})-\rho(T_{exit})}{\rho(T_{bottom})}\cdot w_{exit} \pi R^{2}, \tag{12}\]  
\[w_{exit}\cdot R^{2}=1.8544\times 10^{-3}\;\mathrm{m}\cdot\mathrm{s}^{-1}. \tag{13}\]  
The effective radius of the hydrothermal vent is 0.1 m, and the initial velocity of the vent fluid is 0.2 m s\({}^{-1}\). In this paper, the wall conditions are divided into wall, symmetry and axis. The seafloor is a true wall, so it is set as a wall with a no-slip boundary condition. The remaining computational domain boundaries are set as axial boundary condition and symmetric boundary condition according to geometric characteristics. Therefore, the specific dettings are as follows. At the vent orifice, the velocity-inlet boundary condition is prescribed with constant exit fluid velocity (w\({}_{\mathrm{exit}}\) = 0.2 m/s) and temperature (T\({}_{\mathrm{exit}}\) = 362 \({}^{\circ}\)C). At the bottom of the domain, a no-slip boundary condition is prescribed with a constant temperature (T\({}_{\mathrm{bottom}}\) = 1.74 \({}^{\circ}\)C). At the top of the domain, a pressure-outlet boundary condition is prescribed with a constant temperature (T\({}_{\mathrm{top}}\) = 1.9565 \({}^{\circ}\)C; depending on the background stratification different cases may have different T\({}_{\mathrm{top}}\)'s); here, the fluid is allowed to flow in or out the boundary but with the constraint of maintaining a constant static pressure. The vertical boundary of the domain representing the plume axis is prescribed with an axis boundary condition. The opposing vertical boundary is prescribed with a symmetry boundary condition. The details of the different boundary conditions used in this study are shown in Table 3.  
Since the bottom is a wall with a non-slip boundary condition, the boundary layer naturally exists. The green line in Figure 4 is the corresponding boundary layer. When the realizable k-epsilon model is used for calculation, the standard wall function is chosen as the wall function in this paper. When using the standard wall function, it is necessary to ensure that the first layer of grid nodes is arranged in the logarithmic region. This means that we need to adjust the y+ values (i.e., the vertical distance between the first layer grid nodes and the wall). In the first attempt, we estimated that y+ is equal to 30. After the calculation, the ANSYS Fluent (version 2020 R2) post-processing software is used to determine whether the wall y+ distribution meets the requirements of the model. Otherwise, the grids and calculation need to be redivided.