Assessment of RANS Turbulence Models in Prediction of the Hydrothermal Plume in the Longqi Hydrothermal Field

Abstract

In this paper, the numerical models are selected to simulate the hydrothermal plume based on the water temperature observation data of the Longqi hydrothermal field in the Southwest Indian Ridge (SWIR). Then, the unsteady Reynolds-averaged Navier–Stokes equations are solved to evaluate the performance of the Realizable k-ε (rke) model and the SST k-ω (sst) model in hydrothermal plume simulation. By comparing the calculated results with the Conductivity Temperature Depth (CTD) observation data and the literature results, the difference in prediction performance between the two models is evaluated. Before the numerical simulation, the optimal mesh parameters are determined by considering the grid independence test. The results show that the relative difference of the maximum plume height calculated by the two models is within 5%. Compared with the CTD 05-2, the rke model calculates the root mean square error of the velocity is 0.5081, which is smaller than that of the sst model. In terms of turbulent viscosity, the rke model is in good agreement with reference value in predicting turbulent viscosity. Therefore, the turbulent viscosity distribution calculated by the rke model is more consistent with the plume development process than that calculated by the sst model. In addition, the two models have the same effect on the prediction of turbulent kinetic energy and plume temperature.

1. Introduction

Submarine hydrothermal plume is formed from the dense mixing of high-temperature fluid emitted from submarine hydrothermal vents with surrounding seawater and is important in the global material circulation and heat energy transport [1,2,3,4,5]. At the same time, the dynamic characteristics of hydrothermal plume are the basis of studying the metallogenic law of hydrothermal system, and the study of the dynamic characteristics of hydrothermal plume can promote the exploration of seabed polymetallic sulfide resources [6,7,8].

Table 1. Symbol specification.

Symbol	Explanation	Units
ρbottom	Density at seafloor	kg m−3
Ttop	Temperature at the model top	°C
Tbottom	Temperature at seafloor	°C
Texit	Temperature of vent exit fluid	°C
wexit	Upward velocity of vent exit fluid	m s−1
Qexit	Volume flux issued from vent orifice	m3 s−1
N	Background buoyancy frequency	s−1
Bexit	Buoyancy flux issued from vent orifice	m4 s−3
Zmax	Maximum plume rise height	m
Zneutral	Neutrally buoyant height of the plume	m

Table 1. Symbol specification.

Symbol	Explanation	Units
ρbottom	Density at seafloor	kg m−3
Ttop	Temperature at the model top	°C
Tbottom	Temperature at seafloor	°C
Texit	Temperature of vent exit fluid	°C
wexit	Upward velocity of vent exit fluid	m s−1
Qexit	Volume flux issued from vent orifice	m3 s−1
N	Background buoyancy frequency	s−1
Bexit	Buoyancy flux issued from vent orifice	m4 s−3
Zmax	Maximum plume rise height	m
Zneutral	Neutrally buoyant height of the plume	m

Nowadays, researchers have studied deep-sea hydrothermal plumes in several ways, such as field observations, laboratory experiments, theoretical analytical models, and computational fluid dynamics models. In terms of theoretical models, Morton et al. simplified the dynamic equation of the plume and proposed a classical Morton–Taylor–Turner (MTT) model which calculated the maximum plume rising height of an ideal point source in linear stratification. Then, Morton et al. extended the model to apply to the case of forced plumes with initial momentum [9,10]. Rooney et al. proposed a potential flow integral model to simulate the mixing process of multi-vent plumes. The result showed that the mixing height of the two plumes is lower than the sum of the heights of each isolated plume [11]. In terms of laboratory experiments, Zhang et al. simulated the development process of hydrothermal plumes under laboratory conditions through a linear stratified brine generator and remeasured the coefficients of the scaling law of rising height [12]. He also estimated the maximum turbulent viscosity of the plume. Mirajkar et al. found that the maximum plume rising height was negatively correlated with the environmental buoyancy frequency through plume experiments with different environmental buoyancy frequencies [13]. With the improvement of computing power, the Computational Fluid Dynamics (CFD) models have been widely used in the prediction of hydrothermal plumes due to their excellent repeatability and controllability. Jiang et al. established a numerical model based on the measured data in the hydrothermal field to study the simulation of the hydrothermal plume by a k-ε model. The results showed the variation trend of turbulent kinetic energy, turbulent dissipation rate and turbulent viscosity [14]. Lou et al. simulated the mixing process of two interacting hydrothermal plumes using the k-ε model further. The result showed that maximum turbulent kinetic energy and maximum turbulent dissipation rate of the plumes decrease with the increase in the plume source spacing, while maximum turbulent viscosity changes in the opposite trend [15].

In the current numerical study of hydrothermal plumes, most literature does not discuss the selection of the Reynolds-average Navier–Stokes (RANS) numerical model, but directly uses the model commonly used in fluid analysis to predict hydrothermal plumes. However, we think it is worth studying in selecting a suitable numerical model for hydrothermal plume. If a suitable numerical model can be selected, more accurate prediction results of hydrothermal plume can be obtained. Mokhtarzadeh-Dehghan et al. used different models to simulate the plume of the double jet [16]. The plume motion under different density and background velocity conditions is analyzed and compared with the experimental results of wind tunnel. The results show that the numerical simulation method can effectively simulate the plume. Compared with the calculation results of the standard k-ε model and the RNG k-ε model, it is shown that the model established by the authors can obtain better calculation results than the standard k-ε model. Kumar et al. numerically simulated the plume using the classical k-ε turbulence model and compared the results with Mirajkar’s experimental value of the plume velocity [17]. Then, the authors reported that the root mean square error (RMS) of velocity of the standard k-ε model is smaller than that of the RNG k-ε model, that is, the calculation effect of the standard k-ε model is better. In addition, Mirajkar et al. measured plume velocities in linear stratified environments using flume experiments [18]. However, Kumar et al. did not deeply study the prediction of turbulence statistics by different k-ε models. The prediction results of turbulence statistics and other characteristic parameters of deep-sea hydrothermal plumes by different RANS turbulence models will be evaluated in depth. The comparison of numerical models also exists in other fields of research. Nie et al. evaluated the performance of six low-Reynolds-number k-ε models such as Abid (AB), Lam and Bremhorst (LB), Launder and Sharma (LS), Yang and Shih (YS), Abe Kondoh and Nagano (AKN), and Chang, Hsieh and Chen (CHC) for the prediction of transitional ventilation flow [19]. After comparing the measured values, Nie et al. believed that the LS model had advantages in the prediction of velocity, and the prediction results of the AB model were in the best agreement with the experiment.

In this study, the numerical models are selected to simulate the hydrothermal plume based on the observation data of Longqi hydrothermal field in the Southwest Indian Ridge (SWIR). We use Fluent simulation to study the performance of the rke model and the sst model for forced plume simulation in linear stratified environment. The simulated velocity field, temperature field and rising height are compared with the observation data. Then, the performance of the rke model and the sst model in predicting the dynamic characteristics of deep-sea hydrothermal plume is evaluated. This paper will help to deepen the study of hydrothermal plumes and promote the selection of suitable numerical models for calculation.

2. Materials and Methods

2.1. Turbulence Model

The unsteady Reynolds-averaged Navier–Stokes equations are used to describe fluid motion in a two-dimensional circular hole jet with a single nozzle. The conservation equation is expressed in the form of Cartesian tensor as follows.

Continuity equation:

$$\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho u_i\right)}{\partial x_i}=0.$$

(1)

Momentum equation:

$$\frac{\partial \left(\rho u_i\right)}{\partial t}+\frac{\partial \left(\rho u_i{u}_j\right)}{\partial x_j}=-\frac{\partial p}{\partial x_i}+\frac{\partial {\tau }_{ij}}{\partial x_j}+\frac{\partial }{\partial x_j}\left(-\rho \overline{u_i^{\prime }{u}_j^{\prime }}\right)+\rho g_i.$$

(2)

Energy equation:

$$\frac{\partial \left(\rho E\right)}{\partial t}+\frac{\partial \left[u_j\left(\rho E+p\right)\right]}{\partial x_j}=\frac{\partial }{\partial x_j}\left({\kappa }_{eff}\frac{\partial T}{\partial x_j}\right)+\frac{\partial }{\partial x_j}\left[u_i{\mu }_{eff}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}-\frac{2}{3}{\delta }_{ij}\frac{\partial u_l}{\partial x_l}\right)\right],$$

(3)

where strain rate tensor ${\tau }_{ij}$ can be expressed as

$${\tau }_{ij}=\mu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}-\frac{2}{3}{\delta }_{ij}\frac{\partial u_l}{\partial x_l}\right).$$

(4)

Through the Boussinesq hypothesis, Reynolds stress of Equation (2) can be expressed as

$$-\rho \overline{u_i^{\prime }{u}_j^{\prime }}={\mu }_t\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)-\frac{2}{3}{\delta }_{ij}\left(\rho k+{\mu }_t\frac{\partial u_l}{\partial x_l}\right).$$

(5)

The turbulence model theory is used to close the above equations. At present, the k-ε model and the k-ω model are used to calculate plume. In this paper, the Realizable k-ε (rke) model and the SST k-ω (sst) model are selected for numerical simulation of dynamic characteristics of hydrothermal plume. The explanation of relevant physical quantities of the numerical model of hydrothermal plume in this paper is shown in Table 1.

The turbulent kinetic energy transport equation and turbulent dissipation rate transport equation of the rke model are as follows, described in detail in Appendix A [20]:

$$\frac{\partial \left(\rho k\right)}{\partial t}+\frac{\partial \left(\rho k{u}_j\right)}{\partial x_j}=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right]+G_k+G_b-\rho \epsilon ,$$

(6)

$$\frac{\partial \left(\rho \epsilon \right)}{\partial t}+\frac{\partial \left(\rho \epsilon u_j\right)}{\partial x_j}=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_j}\right]+\rho C_{1}S\epsilon -\frac{\rho C_2{\epsilon }^2}{k+\sqrt{\epsilon v}}+C_{1\epsilon }\frac{\epsilon }{k}{C}_{3\epsilon }{G}_b.$$

(7)

Similarly, the sst model can be written as [21]

$$\frac{\partial \left(\rho k\right)}{\partial t}+\frac{\partial \left(\rho k{u}_j\right)}{\partial x_j}=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right]+P_k-{\beta }^{*}\rho k\omega ,$$

(8)

$$\frac{\partial \left(\rho \omega \right)}{\partial t}+\frac{\partial \left(\rho \omega u_j\right)}{\partial x_j}=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\omega }}\right)\frac{\partial \omega }{\partial x_j}\right]+P_{\omega }-\beta \rho {\omega }^2+2\left(1-F_1\right)\frac{\rho {\sigma }_{\omega 2}}{\omega }\frac{\partial k}{\partial x_j}\frac{\partial \omega }{\partial x_j}.$$

(9)

2.2. Hydrological Background of the Longqi Hydrothermal Field

The SWIR is a typical ultra-slow spreading ridge [22,23]. Figure 1a shows that the SWIR has obvious segmental property, and the SWIR is the best place to study the relationship between magmatic activity, tectonics and hydrothermal circulation of mid-ocean ridges [24]. The Longqi hydrothermal field is located at the junction of the ridge axis and non-conversion discontinuity in the 28th segment of the SWIR, on a hump of the southeast slope of the central axis rift, with a water depth of about 2750 m [25]. Among them, the Longqi-1 hydrothermal field is on the southern side of SWIR segment 28 [6]. It can be seen from Figure 1 that the terrain in the middle and eastern section of the ridge is relatively gentle, the ridge rift is shallow and narrow, the central axis rift in the western section of the ridge is wide and deep, and the water depth changes greatly. The topographic structure with fissures provides channels for seawater infiltration into the deep crust and upwelling of hydrothermal fluids and provides favorable topographic environment for hydrothermal activities.

The first active hydrothermal field was discovered in the 28th ridge segment of SWIR during the 19th Chinese Oceanic voyage in 2007, and it was named the Dragon Flag Hydrothermal Field (DFHF, 49.65° E, 37.78° S) [25,26,27]. The Longqi hydrothermal field has three hydrothermal vent zones: the N zone, the M zone and the S zone. During the 35 Ocean voyages of China from 2014 to 2015, Jiaolong was used to make detailed near-bottom observations of the hydrothermal field of Longqi, and it was determined that the Longqi hydrothermal field is one of the large-scale hydrothermal vents of the mid-ocean ridge [6,22,28]. Two active hydrothermal zones about 500 m apart, called the S and the M zones, have been found at depths of 2750–2780 m in the Longqi-1 hydrothermal field [25]. Observations made by the researchers in 2016 using a HOV (Human-Occupied Vehicle) indicated that there are six high-temperature vents in the M zone and three in the S zone. At the same time, five large diffuse areas were found in the S zone [28]. Tao et al. investigated the hydrothermal zone of Longqi-1 (49.645° E~49.650° E, 37.780° S~37.785° S) and found three active spouts (DFF3, DFF5 and DFF20) in the S zone (Figure 1b and

Table 2. Active hydrothermal vents of Longqi-1 (S zone). Specific data come from Tao et al. [6].

Vent	Dive-Sampler	Max Texit (°C)	Lon (° E)	Lat (° S)
DFF 3	JL 89-CGT-B/C	352	49.649503	37.782617
DFF 5	-	146	49.649910	37.783365
DFF 20	JL 96-CGT-D/E	362	49.649093	37.783563

) [6,28,29]. It was recently found that a melt zone in the lithospheric mantle ~13 ± 2 km below the sea floor is the most likely source of heat for the Longqi-1 hydrothermal field, and the geochemical composition of the Longqi-1 hydrothermal fluids implied that the path of the Longqi-1 hydrothermal circulation may be longer than any other hydrothermal system studied to date [6]. Figure 1c shows the distribution of CTD stations marked by our team on the fourth leg of the 30th voyage. Among them, the hydrothermal plume near the CTD 05-2 station is the research object of this paper.

Miniature Automatic Plume Record (MAPR) is a tool used in ocean observations to find hydrothermal vents. MAPR is equipped with CTD to carry out launching operations with a maximum diving depth of 2835 m, entry time of 16:25, end time of 18:46, and exit time of 21:34, for a total underwater operation of 5 h and 09 min. The results are shown in

Table A1. Summary of MAPR operation in Section IV of Ocean 30 Voyage (CTD05-2-MAPR).

Date/Time (dd/mm/yyyy hh:mm:ss)	Temp (°C)	Depth (m)	Neph (V)
22/04/2014 16:39:01	19.77275	95.71033	0.18492
22/04/2014 16:39:06	19.68329	99.31679	0.17333
22/04/2014 16:39:11	19.57919	103.26665	0.24154
22/04/2014 16:39:16	19.54152	104.64050	0.17740
22/04/2014 16:39:21	19.47655	104.29704	0.17532
...	...	...	...
22/04/2014 18:07:06	2.81011	2003.16878	0.16054
22/04/2014 18:07:11	2.80706	2005.89159	0.16279
22/04/2014 18:07:16	2.80629	2007.93368	0.15821
22/04/2014 18:07:21	2.80400	2010.14591	0.16270
22/04/2014 18:07:26	2.80324	2013.54931	0.15950
...	...	...	...
22/04/2014 19:30:01	2.28663	2626.80683	0.16581
22/04/2014 19:30:06	2.28588	2627.65524	0.19659
22/04/2014 19:30:11	2.28588	2625.27968	0.16512
22/04/2014 19:30:16	2.28512	2624.60095	0.17541
22/04/2014 19:30:21	2.28437	2622.73442	0.16702

and Figure A1 of Appendix B. In addition,

Table A2. Summary table of CTD stations in the Longqi hydrothermal field.

Voyage Segment	Station No.	Lon (° E)	Lat (° S)	Depth (m)
30IV	30IV-SWIR-S024-CTD02-1	49.646165	37.782502	2807
30IV	30IV-SWIR-S024-CTD03-1	49.649092	37.783305	2812
30IV	30IV-SWIR-S024-CTD03-3	49.649697	37.782720	2804
30IV	30IV-SWIR-S024-CTD04	49.648712	37.782807	2808
30IV	30IV-SWIR-S024-CTD05-1	49.651808	37.782900	2838
30IV	30IV-SWIR-S024-CTD05-2	49.649380	37.783912	2818
30IV	30IV-SWIR-S024-CTD06	49.650255	37.781202	2817

and

Table A3. CTD05-2 Station information collected.

Depth (m)	Temp (°C)	Sal (PSU)	Depth (m)	Temp (°C)	Sal (PSU)
2000.497	2.4754	34.8941	2729.825	1.9602	34.9823
2001.278	2.4759	34.8939	2729.537	1.9603	34.9821
2001.482	2.4760	34.8936	2728.868	1.9605	34.9820
2001.447	2.4757	34.8936	2727.958	1.9603	34.9821
2001.291	2.4780	34.8930	2726.954	1.9601	34.9821
2001.055	2.4794	34.8926	2726.075	1.9599	34.9821
2001.022	2.4790	34.8927	2725.580	1.9598	34.9822
2001.478	2.4779	34.8932	2725.594	1.9597	34.9823
2002.352	2.4779	34.8931	2725.886	1.9595	34.9825
2003.383	2.4725	34.8947	2726.121	1.9596	34.9822
2004.332	2.4731	34.8946	2725.915	1.9598	34.9821
2005.050	2.4769	34.8932	2725.224	1.9598	34.9821
2005.462	2.4743	34.8941	2724.337	1.9599	34.9820
2005.580	2.4739	34.8942	2723.440	1.9596	34.9821
2005.476	2.4740	34.8942	2722.551	1.9594	34.9821
2005.369	2.4756	34.8937	2721.560	1.9594	34.9820
2005.603	2.4741	34.8942	2720.476	1.9595	34.9818
2006.354	2.4735	34.8944	2719.486	1.9595	34.9819
2007.433	2.4738	34.8943	2718.892	1.9593	34.9821
2008.638	2.4739	34.8944	2719.001	1.9591	34.9823

in Appendix C list the CTD stations and their observations.

The temperature and turbidity data collected from Conductivity Temperature Depth (CTD) stations are processed and analyzed [30,31]. The results show that there are obvious hydrothermal activities near CTD 02-1, CTD 03-3 and CTD 05-2 observation stations and the overall turbidity anomaly variation range is within 0.05–0.15 NTU (Figure 2a,b). The water depth of the CTD 05-2 station is about 2730 m. The yellow area (2420~2540 m) has obvious turbidity anomaly, and it is speculated that this area (190~310 m from the bottom) is affected by the neutral buoyant layer of the hydrothermal plume (Figure 3a). According to the corresponding temperature value, the seawater temperature gradient in the sea area around the observation station is determined to be 0.0002165 K m⁻¹.

The thermohaline environment and spatial scale of the hydrothermal plume in the process of movement are mainly affected by temperature. Based on the water depth (2730 m), pressure (about 27 MPa) and salinity (34.98 PSU) information of the CTD 05-2 station in the Longqi hydrothermal field, seawater temperature is taken as a functional independent variable, and density and specific heat capacity of seawater are calculated using the calculation method of seawater physical property parameters provided by other scholars [32]. Figure 3b shows the variation of density and specific heat capacity of seawater with temperature, respectively. The image shows that the influence of seawater temperature on density is greater than that on specific heat capacity. In the bottom seawater environment (temperature below 4 °C), the density curve is approximately a straight line, and the rate of change in seawater density to temperature is 0.0786 kg m⁻³ K⁻¹.

2.3. 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},$$

(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 \left(T_{bottom}\right)}\cdot \frac{\partial \rho \left(z\right)}{\partial z}}=\sqrt{\frac{g}{\rho \left(T_{bottom}\right)}\cdot \frac{\partial \rho }{\partial T}\cdot \frac{\partial T}{\partial z}},$$

(11)

$$B_{exit}=g\frac{\rho \left(T_{bottom}\right)-\rho \left(T_{exit}\right)}{\rho \left(T_{bottom}\right)}\cdot w_{exit}\pi R^2,$$

(12)

$$w_{exit}\cdot R^2=1.8544\times {10}^{-3}\mathrm{m}\cdot {\mathrm{s}}^{-1}.$$

(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⁻¹.

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_exit = 0.2 m/s) and temperature (T_exit = 362 °C). At the bottom of the domain, a no-slip boundary condition is prescribed with a constant temperature (T_bottom = 1.74 °C). At the top of the domain, a pressure–outlet boundary condition is prescribed with a constant temperature (T_top = 1.9565 °C; depending on the background stratification different cases may have different T_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. Different boundary conditions used in the study.

Type	Vent Orifice	Bottom (Wall)	Top	Sides
Boundary Condition	Velocity Inlet	No-Slip Boundary	Pressure Outlet	Symmetry Boundary
Velocity	wexit = 0.2 m/s	No-Slip (wexit = 0)	Flow in or out the boundary	Zero Gradient *
Temperature	Texit = 362 °C	Tbottom = 1.74 °C	Ttop = 1.9565 °C	Zero Gradient *
Dynamic Pressure	p = 0	p = 0	Constant static pressure	Zero Gradient *

* Zero gradient means that the normal gradient of the physical quantity at the boundary is zero.

.

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. The y+ values of this paper are in the range of 30~100, indicating that the boundary layer grid nodes are reasonably arranged.

2.3.2. Grid Division and Independence Test

In the actual environment, there are interactions between hydrothermal plumes. To effectively simulate the plume formed by hydrothermal vent groups, a simplified scheme of equivalent vents can be adopted [33]. Specific operations: (1) We assume an axisymmetric three-dimensional space and transform it into a two-dimensional plane computation domain with axisymmetric nozzle. (2) According to the flow equivalence, multiple adjacent active spouts in the active spout group are equivalent to an ideal circular spout. Figure 5a,b shows the detailed mesh size and its partitioning: horizontal radial radius is 1000 m, vertical height is 1000 m, effective radius of vent orifice is 0.1 m, and vent orifice is 2 m above the sea floor.

The vent orifice is the initial places where hydrothermal fluids erupt, and in order to effectively capture the fine flow of the plume, the mesh of the vent orifice should be treated. Meanwhile, Suzuki et al. [34] also showed that turbulent mixing near the vent orifice can be correctly reproduced when the mesh size of the area near the vent orifice should be less than D/10 (D is the diameter of the vent orifice). We use the Triangles tool in Fluent to draw the grid. First, the vent orifice radius is evenly divided into 10 equal parts, and then 1/10 of the vent orifice radius is the smallest mesh size (the smallest mesh size is 0.01 m). The grid is set to extend from the vent orifice to the computational domain boundary at a constant expansion rate of 1.02 (Figure 5a).

In this paper, grid independence test is carried out to verify rationality of grid di-vision, reduce the errors caused by poor grid quality in numerical simulation, and effectively simulate the plume movement in the vent orifice. The mesh of vent orifice radius is divided into seven different scales, namely, 3, 5, 7, 9, 11, 13 and 15 (Figure 6 and

Table 4. Grid independence test.

Mesh Index	Mesh Partition	Zneutral (m)	Zmax (m)	Zneutral/Zmax	Total Number of Elements
1	3	243.427765	333.467255	0.729990	43,968
2	5	242.232590	331.873627	0.729894	49,465
3	7	239.443848	326.295929	0.733824	53,378
4	9	239.045456	326.295929	0.732603	56,520
5	11	242.630981	325.897522	0.744501	59,026
7	13	240.016435	325.365413	0.737683	61,832
8	15	239.431358	325.464651	0.735660	63,484

).

The maximum plume rising height is the vertical height when the vertical velocity of the plume centerline becomes 0 for the first time. The neutrally buoyant plume height is the vertical height when the density of the plume centerline is equal to the density of the background seawater for the first time.

It is obvious that the results of Scheme 1 and Scheme 2 are significantly different from those of the other five schemes, which may be caused by insufficient grid detail. The results of the last five schemes of calculation schemes are generally similar, and the differences between the last five schemes are significantly smaller than the differences between the first two schemes. For the maximum plume rising height, the relative difference of the calculated results is less than 0.5%. For the neutrally buoyant plume height, the relative difference among all groups is 1% except Scheme 5. This paper considers that there is a possibility of calculation error in Scheme 5. As can be seen from Figure 6 and Table 4, the corresponding calculation results tend to be stable as the degree of grid division increases. In addition, the increase in the total number of grids in the computing domain also tends to be stable. Therefore, to save time and resources, the design of Scheme 3 (vent orifice 7 equal division) is selected for grid division. Finally, the entire computing domain is discretized by 53,378 triangular elements, as shown in Figure 5b.

2.3.3. Numerical Simulation of Working Conditions

Buoyancy frequency and buoyancy flux are the key factors affecting the movement of hydrothermal plumes. Taking the environmental reference condition (case2) as the control group, environmental factors such as the initial velocity of the vent orifice fluid, the temperature of the vent orifice fluid and the temperature gradient of the environment are changed to observe the influence of those on the dynamic characteristics of hydrothermal plume. To evaluate the ability of different models to simulate hydrothermal plumes, the rke model and the sst model are selected for comparative analysis, and 14 groups of numerical simulation conditions are set. See

Table 5. Working conditions group of numerical simulation.

Simulation Run (rke and sst)	Texit (°C)	Ttop (°C)	wexit (m s−1)	Bexit (m4 s−3)	N (s−1)
case 1	362	1.957	0.1	0.008183	0.000401
case 2 *	362	1.957	0.2	0.016365	0.000401
case 3	362	1.957	0.4	0.032729	0.000401
case 4	362	2.500	0.2	0.016365	0.000751
case 5	362	3.000	0.2	0.016365	0.000967
case 6	250	1.957	0.2	0.008248	0.000401
case 7	300	1.957	0.2	0.010119	0.000401

* Reference environment group.

for specific simulation parameter settings. In addition, temperature gradient initialization is used in the simulation of each set of conditions.

In order to investigate the process of hydrothermal plume from eruption to quasi-steady state, the calculation time of the simulation is set as follows. The total simulation time of each condition is 4t*, where t* is the buoyancy time scale, defined as t* = 2π/N. We use a second-order implicit time integration scheme in our simulations. An implicit scheme is unconditionally stable with respect to time step size. At the same time, a time step of 0.5 s is selected so that the calculation in a single time step can reach convergence within 20 iterations (that is, the total residual of each physical quantity in two successive iterations is less than 10 × 10⁻³, and the total residual of energy is less than 10 × 10 ⁻⁶). In addition, because the simulated physical time is short relative to the rotation period of the Earth, the effect of the Coriolis force is ignored. When the simulation time is 3t*, the curves of each variable tend to be stable. Therefore, readers can judge the convergence speed and stability (that is, the asymptotic stability of numerical calculations) by observing the convergence of the iterative residual in each working condition.

In this paper, the computational fluid dynamics software ANSYS Fluent (version 2020 R2) based on the finite volume method was used to numerically solve the above turbulence models and the specific settings for hydrothermal plume modeling. In order to reduce the numerical dissipation and ensure the rationality of the numerical model results, the following settings were employed in this paper. The highly accurate third-order MUSCL (monotone upstream-centered schemes for conservation laws) scheme was used for spatial discretization for all the convection terms in the momentum, energy and turbulence equations. The second-order accurate central-differencing scheme was used for spatial discretization for the diffusion terms. We chose the PRESTO! (Pressure Staggering Option) scheme as the discretization method of pressure. The PISO (pressure implicit with splitting of operators) scheme was adopted in the pressure–velocity coupling scheme. Second-order implicit scheme was used for temporal discretization.

3. Results and Discussion

3.1. Temperature Field and Velocity Field

3.1.1. Temperature Field (T)

In the numerical simulation, the temperature gradient is used for initialization calculation. Figure 7 shows the temperature distribution simulated by different turbulence models in a quasi-steady state. When the hydrothermal plume receives vertical momentum, hot fluids mix violently with cold seawater. As hydrothermal plume gradually moves away from the vent orifice and entrails the cooler fluid, its buoyancy effect diminishes gradually and its temperature drops. After the neutrally buoyant plume height, the temperature contour line flattens out, and the temperature of hydrothermal plume gradually approaches the ambient temperature. There are some differences in the temperature distribution between the two models. At the seafloor, the variation interval of the temperature contour of the rke model is smaller than that of the sst model. At the maximum plume rising height, there are some fluctuations in the temperature contour of the rke model, while the temperature contour of the sst model is too idealized.

The relationship between temperature, depth and time during detection is shown in Figure A1 in Appendix B. The relationship between temperature and depth is obtained. The plume temperature calculated by the rke model and the sst model is analyzed based on this temperature. When z/Z_max is less than 0.1, the temperature of the hydrothermal plume decreases rapidly from 362 °C to 4 °C (Figure 8). The variation trend of the temperature curve shows that the simulated temperature variation of the plume is consistent with the law of plume movement. According to the simulation results of environmental conditions (case2), the RMS of temperature of the rke model is 0.0408 and that of the sst model is 0.0406. The root mean square error of temperature and the variation trend of the curve show that the numerical simulation of hydrothermal plume in this paper is in agreement with reality.

In the same working conditions, compared with the reference value, the calculation temperature of the rke model is higher than that of the sst model. But there is relatively little difference between the two models. In conclusion, the prediction effects of the rke model and the sst model are similar in terms of hydrothermal plume temperature.

3.1.2. Velocity Field (w)

When the hydrothermal plume receives vertical momentum, its velocity increases rapidly from the initial velocity to the maximum velocity. Due to the fluid’s gravity, viscous effect, and entrainment, the buoyancy effect diminishes gradually. Its velocity begins to decrease until it becomes zero at the maximum plume rising height. Figure 9 shows the velocity distribution simulated by different turbulence models in a quasi-steady state. The velocity vectors of the two models are obviously different in the near-wall region. Compared with the rke model, the sst model carries out too much calculation on the boundary layer of the wall. The upwelling process of hydrothermal plumes is mainly turbulent entrainment, mixing and diffusion, but it is unlikely to have too much contact with the seafloor boundary layer.

Since the plume velocity was not measured in this observation, the plume velocity simulated by Lou et al. in the Daxi hydrothermal field of the Carlsberg Ridge in a similar environment was selected as the comparison standard [35]. On this basis, the plume velocity calculated by the rke and sst models is analyzed. Figure 10 shows that the vertical velocity of the hydrothermal plume increases from the initial velocity of 0.2 m s⁻¹ to 1 m s⁻¹ or even 2 m s⁻¹ in the area near the vent. This trend is consistent with the movement of hydrothermal plumes. According to the simulation results of environmental conditions (case2), the RMS of velocity of the rke model is 0.5081 and that of the sst model is 0.5122. The comparison of data again shows that the numerical simulation of hydrothermal plume in this paper is reasonable.

In the same working conditions, compared with the reference value, the calculation velocity of the rke model is larger than that of the sst model, and the rke model has a smaller RMS error than the sst model. In addition, the sst model has a transitional prediction of the seabed velocity. To sum up, the rke model is better at predicting the velocity of hydrothermal plume.

3.2. Turbulence Statistics

3.2.1. Turbulent Viscosity (μ_t, ν_t)

As can be seen from Figure 11, the plume turbulent viscosity presents a mushrooming spatial distribution, and the value decays outward from the center of the plume region and the axis of the plume center, and the strongest turbulence mixing exists in the plume region. The maximum turbulent viscosity calculated by rke_case2 and sst_case2 are 100.4 Pa s and 95.4 Pa s. It can be seen from the model itself that the rke model strengthens the ability of standard k-ε model to predict the flow field of rotating flow, and the sst model strengthens the prediction ability of the standard k-ω model in the boundary layer. However, according to the definition, the sst model adopts the algorithm of the standard k-ε model to calculate the jet region in the range of plume stem region. The hydrothermal plume belongs to local axisymmetric flow [36], and there are vortices, enrolling and rotating flow phenomena. Therefore, the turbulent viscosities calculated by the two turbulence models are not the same.

The changes in turbulent viscosity of the vent orifice center line in different working conditions along with the height of hydrothermal plume rising are shown in Figure 12a. The trends calculated by the two models are not the same. When the height z is less than 0.8 Z_max, the turbulent viscosity calculated by the rke model is larger than that calculated by the sst model. Compared with the sst model, the calculated results of the rke model are closer to those of Jiang et al. [14].

Through dimensional analysis, it is concluded that there is a certain functional relationship between turbulent viscosity, buoyancy frequency and buoyancy flux. Figure 12b shows the linear fitting of the maximum turbulent viscosity of the vent orifice center line with the buoyancy flux and buoyancy frequency. The turbulent viscosities scaling law of the rke model and the sst model are as follows:

$$\{\begin{array}{l}{v}_{t,\max ,rke}=0.016B_{exit}^{1/2}\cdot N^{-1/2}\\ v_{t,\max ,sst}=0.014B_{exit}^{1/2}\cdot N^{-1/2}\end{array}.$$

(14)

Zhang et al. obtained the scaling law of turbulent viscosity through plume experiments, and its scaling coefficient was 0.030, an order of magnitude consistent with that in this paper [12]. We believe that the turbulent viscosity scaling law of the hydrothermal plume is reasonable considering the difference between the experimental environment and the simulated environment. From the numerical value of turbulent viscosity, compared with the sst model, the turbulent viscosity calculated by the rke model is larger. There is intense enrolling and vorticity phenomena in the hydrothermal plume, and there are large flows in and below the neutral buoyant layer, which generates large internal friction force and increases the turbulent viscosity. Therefore, the rke model can predict the turbulent viscosity of hydrothermal plume more reasonably.

3.2.2. Turbulent Kinetic Energy (k)

As can be seen from Figure 13, the turbulent kinetic energy of the hydrothermal plume presents a mushroom-like spatial distribution. Its value reaches its maximum in the plume stem area near the vent orifice and decreases sharply along the radial direction. The maximum turbulent kinetic energy calculated by rke_case2 and sst_case2 are 0.069 J kg⁻¹ and 0.083 J kg⁻¹. According to the computational relationship between turbulent kinetic energy and velocity, there is a positive correlation between the velocity field and turbulent kinetic energy. When the hydrothermal plume diffuses laterally in the neutrally buoyant plume height, there are strong anisotropic flows such as vortices. The velocity contour line in Figure 13 indicates that the rke model is sensitive to velocity variation in the far field region, while the sst model is not suitable for free shear flow in the late plume period because of its large magnitude variation in velocity prediction. Compared with the previous turbulent kinetic energy nebulae distribution, the prediction error of the sst model on the transverse diffusion of the plume neutral buoyancy layer is larger [12,14,15].

The changes in turbulent kinetic energy of nozzle center line with plume height under different working conditions are shown in Figure 14a. Since there is no measurement data of turbulent kinetic energy, the conclusion of Zhang et al. was taken as reference [12]. The curves in the area near the vent orifice show that the difference between turbulent kinetic energy calculated by the rke model and the sst model is small, and the change trends of the two models are similar.

Similarly, it is determined that there is a functional relationship among turbulent kinetic energy, buoyancy frequency and buoyancy flux by dimensional analysis. Figure 14b shows the linear fitting of maximum turbulent kinetic energy of the vent orifice center line with the buoyancy flux and buoyancy frequency. The turbulent kinetic energy scaling law of the rke model and the sst model is as follows:

$$\{\begin{array}{l}{k}_{\max ,rke}=22.646B_{exit}^{1/2}\cdot N^{1/2}\\ k_{\max ,sst}=24.291B_{exit}^{1/2}\cdot N^{1/2}\end{array}.$$

(15)

From the numerical analysis of maximum turbulent kinetic energy, under the same working conditions, the maximum turbulent kinetic energy calculated by the sst model is slightly greater than that calculated by the rke model, but the difference is not significant. Except for the difference in turbulent kinetic energy at transverse diffusion, the numerical values of turbulent kinetic energy predicted by the rke model and the sst model are similar to the curve change trend, indicating that these two models have similar effects on the prediction of turbulent kinetic energy of hydrothermal plume.

3.2.3. Turbulent Dissipation Rate (ε, ω)

Turbulent dissipation rate $\epsilon$ and specific turbulent dissipation rate $\omega$ represent the dissipation rate of the rke model and the sst model. The relationship between the two is as follows:

$$\omega =\frac{\epsilon }{C_{\mu }\cdot k}.$$

(16)

As can be seen from Figure 15, the turbulent kinetic energy dissipation rate of the hydrothermal plume presents a mushroom-like spatial distribution. Its value reaches its maximum in the plume stem area near the vent orifice and decreases sharply along the radial direction. Compared with the turbulent kinetic energy, the radial attenuation of the turbulent dissipation rate has a greater degree of order of magnitude change. After the conversion of the two turbulent dissipation rates, the maximum turbulent dissipation rates calculated by rke_case2 and sst_case2 are 0.2 W kg⁻¹ and 0.18 W kg⁻¹. Further comparing with the results of Zhang et al., the dissipation rate distribution of the sst model is too exaggerated [12]. Since turbulent dissipation rate represents the rate at which turbulent kinetic energy is converted into heat energy under the action of molecular viscosity, the outline of turbulent kinetic energy and turbulent dissipation rate cloud map is generally similar. The rke model improves the calculation accuracy of the rotating flow by improving the structure of ε transport equation. The prediction of rotational flow and free shear flow in the far field region by the sst model is not as accurate as that by the rke model, which leads to the high deviation of the sst model in predicting turbulent dissipation rate of the neutral buoyant layer.

3.3. Feature Height and Entrainment Coefficient

3.3.1. Feature Height (Z_max, Z_neutral)

Based on the vertical velocity and inner and outer density difference in the hydrothermal plume, the maximum plume rising height of the plume and the neutrally buoyant plume height are obtained (

Table 6. Hydrothermal plume simulation results of rke model and sst model.

Simulation Run		Bexit (m4 s−3)	Basymp (m4 s−3)	Zneutral (m)	Zmax (m)	Zmax,MTT (m)	νt,max (m2 s−1)	kmax (m2 s−2)	εmax (m2 s−2) or ω (s−1)
rke	case1	0.008183	0.000871	215.931	298.798	228.012	0.088830	0.045862	0.098373
	case2 *	0.016365	0.001742	237.444	324.296	271.154	0.096543	0.064709	0.182601
	case3	0.032729	0.003484	285.216	386.447	322.458	0.098900	0.101572	0.292038
	case4	0.016365	0.001742	154.578	213.540	169.320	0.093237	0.073890	0.213846
	case5	0.016365	0.001742	121.113	164.538	140.079	0.083291	0.064936	0.183805
	case6	0.008248	0.001200	216.910	299.595	247.052	0.091085	0.045349	0.088810
	case7	0.010119	0.001442	227.484	315.133	258.649	0.095439	0.053072	0.120054
sst	case1	0.008182	0.000871	204.377	280.073	228.012	0.069245	0.048405	29.112717
	case2 *	0.016365	0.001742	237.045	320.710	271.154	0.090891	0.075377	38.509010
	case3	0.032729	0.003484	288.355	391.499	322.458	0.100053	0.104194	47.898373
	case4	0.016365	0.001742	149.798	209.158	169.320	0.067635	0.075532	38.452026
	case5	0.016365	0.001742	124.301	170.514	140.079	0.060755	0.075374	38.525658
	case6	0.008248	0.001200	231.468	313.539	247.052	0.075987	0.046407	29.675610
	case7	0.010119	0.001442	240.631	326.288	258.649	0.082737	0.055771	32.198055

* Reference environment group.

). The plume neutral buoyancy zone inferred from the detection data indicates that the neutrally buoyant plume heights calculated by rke_case2 and sst_case2 are consistent with the observation results of this scientific investigation. The rke model and the sst model have a height ratio Z_neutral/Z_max of 0.72~0.74 in all working conditions, with an average value of 0.7304. These are close to the empirical value of 0.761 of the classical plume theoretical model [9]. The results show that RANS numerical simulation method can effectively reproduce the macroscopic dynamic process of hydrothermal plume in the Longqi hydrothermal field.

As there is a certain gap between the numerical simulation and the actual environment buoyancy flux, the calculation method of the initial buoyancy flux is modified [37]. Buoyancy flux $B_0$ can be calculated in two ways: the initial buoyancy flux $B_{exit}$ and the corrected progressive buoyancy flux $B_{asymp}$. As can be seen from Equation (18), when the temperature $T_{bottom}$ is 1.74 °C, the seabed density ${\rho }_{bottom}$ is 1038 kg m⁻³ and the thermal expansion coefficient ${\beta }_{bottom}$ is 8 × 10⁻⁵ K⁻¹.

$$B_{asymp}=g{\beta }_{bottom}\left(T_{exit}-T_{bottom}\right)Q_{exit},$$

(17)

$${\beta }_{bottom}=\frac{1}{V}{\left(\frac{\partial V}{\partial T}\right)}_{P,S}=-\frac{1}{\rho }{\left(\frac{\partial \rho }{\partial T}\right)}_{P,S}.$$

(18)

The empirical coefficient of the plume theoretical model was determined by many experiments, and its value was generally recommended to be 3.76 [9]. Figure 16a shows the linear fitting of the maximum plume rising height with the buoyancy flux and buoyancy frequency. Figure 16b shows the linear fitting of the neutrally buoyant plume height, the buoyancy flux and the buoyancy frequency of the hydro-thermal plume.

Table 7. Scaling law of maximum plume rising height and neutrally buoyant plume height.

Scaling Law of Zmax	Scaling Law of Zneutral
Zmax,rke = 2.68 (Bexit N−3)1/4	Zneutral,rke = 1.95 (Bexit N−3)1/4
Zmax,rke = 4.59 (Basymp N−3)1/4	Zneutral,rke = 3.34 (Basymp N−3)1/4
Zmax,sst = 2.69 (Bexit N−3)1/4	Zneutral,sst = 1.97 (Bexit N−3)1/4
Zmax,sst = 4.61 (Basymp N−3)1/4	Zneutral,sst = 3.39 (Basymp N−3)1/4

shows the scaling law of feature height obtained according to the working conditions of the rke model and the sst model.

The MTT theoretical model predicts the result based on an ideal point source, while the numerical simulation predicts the result based on a certain vent radius. The results obtained by using the progressive buoyancy flux are closer to the results of the MTT theoretical model than those obtained by using the initial buoyancy flux. Zhang et al. proposed that the empirical coefficient of the maximum plume rising height was 3.56 through the laboratory experiment of plume, which was 5.3% higher than the empirical value of the theoretical model [12]. When the initial buoyancy flux is adopted, the scaling law coefficients of neutrally buoyant plume height and maximum plume rising height of rke model and sst model are like the results of Lou et al. [35]. Their fitting stability R² is greater than 0.99. When the progressive buoyancy flux is used, the deviation of the fitting coefficient between the rke model and theoretical model is smaller than that between the sst model and the theoretical model. The heights calculated by the rke model and the sst model are not much different. The rke model and the sst model have equivalent prediction effects on the characteristic height of hydrothermal plume, and they can reproduce the macroscopic kinetic process of hydrothermal plume well.

3.3.2. Entrainment Coefficient (α_e)

The entrainment strength of hydrothermal plume is measured by entrainment coefficient. According to analyzing the calculation method of the entrainment coefficient, the calculation formula is defined as [34,38,39]

$${\alpha }_e(z)=\frac{dM(z)}{dz}\cdot \frac{1}{2\pi {\rho }_0\cdot w\cdot r}.$$

(19)

Figure 17a shows that the entrainment coefficients of the rke model and the sst model have a similar variation trend. When z/Z_max is greater than 0.3, the entrainment coefficient increases rapidly from 0 to a larger value. When z/Z_max is 0.03 to 0.60, the entrainment coefficient is stable within 0.10~0.11. When z/Z_max is 0.60 to 0.75, the entrainment coefficient decreases rapidly from a stable value to 0. When z/Z_max is less than 0.5, the average entrainment coefficients calculated by the rke model and the sst model are analyzed. The result of the rke model is 0.10760, which is slightly higher than that of the sst model (the sst model result is 0.10695). Figure 17b shows the difference between the entrainment coefficient calculated by the rke model and the sst model under environmental reference conditions. We take the average entrainment coefficient α_e is 0.104 of Lou et al. as a reference and compare the entrainment coefficient calculated by the rke model and the sst model [35]. The calculation error of rke_case2 is 2.71%, and that of sst_case2 is 1.94%. In terms of the prediction of the entrainment coefficient of hydrothermal plume, the resulting curve of the sst model tends to be gentle, while the velocity change and enrolling effect of the hydrothermal plume in the plume stem area are relatively intense.

Lou et al. introduced the Richardson number (Ri) to calculate the entrainment coefficient [35]. However, when z/Z_max is greater than 0.6, there is a slight deviation in the results of this calculation method. To better evaluate the accuracy of the entrainment coefficients of the rke model and the sst model, the Froude number (Fr) is introduced to improve the calculation method of entrainment coefficient [40].

$$Fr={\left[\frac{w_0^2\cdot {\rho }_a}{gd({\rho }_a-{\rho }_0)}\right]}^{1/2},$$

(20)

$${\alpha }_{Fr}={\alpha }_j-\left({\alpha }_j-{\alpha }_p\right){\left(\frac{F{r}_p}{Fr}\right)}^2.$$

(21)

List et al. determined that the entrainment coefficient of pure jet α_j is 0.0742 and the entrainment coefficient of pure plume α_p is 0.1178 [41]. According to the above equation, we calculate Froude number Fr_p is 4.71879. Figure 17c,d show good self-consistency between α_Fr and α_e. It shows that the entrainment coefficients obtained by the two methods in this paper are correct. The normal range of entrainment coefficient is 0.10~0.16 [42,43,44,45]. The results show that these two models can effectively reproduce the enrolling process of hydrothermal plume.

4. Conclusions

This paper evaluates the numerical results of the rke model and the sst model on the dynamic characteristics of deep-sea hydrothermal plume in the Longqi hydrothermal field. The calculated results are compared with the CTD 05-2 observation data, Kumar’s results and Jiang’s results. Therefore, the following conclusions can be drawn.

The turbulent viscosity obtained by the rke model is larger than that obtained by the sst model, and its trend is closer to the reference. These results show that the rke model can predict the turbulent viscosity of hydrothermal plume more reasonably. Except for the difference in the turbulent kinetic energy at the horizontal diffusion, the maximum turbulent kinetic energy calculated by the two models is not much different, and the corresponding curves have a similar change trend. These results show that the rke model and the sst model have equivalent effect on the prediction of turbulent kinetic energy of hydrothermal plume. The high deviation of the sst model in predicting turbulent dissipation rate of neutral buoyant layer and its prediction results are not as accurate as those of the rke model.

In the vertical velocity of the center line of the hydrothermal plume, the velocity of the RMS of the rke model is smaller than that of the sst model. This shows that the rke model can better predict the velocity of hydrothermal plume. Compared with the data of the CTD 05-2 observation station, the temperature of the RMS of the rke model is slightly larger than that of the sst model. But overall, both models calculate temperatures that are close to the observed data. Therefore, the rke model and the sst model can predict the hydro-thermal plume temperature in a similar way.

The results of characteristic heights show that both the rke model and the sst model can predict macroscopical processes of hydrothermal plumes well. Overall, the two commonly used RANS models accurately capture the macroscopic physical characteristics of buoyancy plumes in stratified environments. But the rke model is more suitable for numerical study of hydrothermal plume than the sst model. At present, due to the lack of experimental data, it is impossible to compare the entrainment coefficient quantitatively. Observation and experiment should be required in the Longqi hydrothermal field in the future.