Physical controls on mixing and transport within rising submarine hydrothermal plumes: A numerical simulation study

Deep-Sea Research I 92 (2014) 41–55

Contents lists available at ScienceDirect

Deep-Sea Research I journal homepage: www.elsevier.com/locate/dsri

Physical controls on mixing and transport within rising submarine hydrothermal plumes: A numerical simulation study Houshuo Jiang n, John A. Breier Department of Applied Ocean Physics and Engineering, Woods Hole Oceanographic Institution, Woods Hole, MA 02543, USA

art ic l e i nf o

a b s t r a c t

Article history: Received 30 December 2013 Received in revised form 28 April 2014 Accepted 14 June 2014 Available online 27 June 2014

A computational fluid dynamics (CFD) model was developed to simulate the turbulent flow and species transport of deep-sea high temperature hydrothermal plumes. The model solves numerically the density weighted unsteady Reynolds-averaged Navier–Stokes equations and energy equation and the species transport equation. Turbulent entrainment and mixing is modeled by a k–ε turbulence closure model. The CFD model explicitly considers realistic vent chimney geometry, vent exit fluid temperature and velocity, and background stratification. The model uses field measurements as model inputs and has been validated by field data. These measurements and data, including vent temperature and plume physical structure, were made in the ABE hydrothermal field of the Eastern Lau Spreading Center. A parametric sensitivity study based on this CFD model was conducted to determine the relative importance of vent exit velocity, background stratification, and chimney height on the mixing of vent fluid and seawater. The CFD model was also used to derive several important scalings that are relevant to understanding plume impact on the ocean. These scalings include maximum plume rise height, neutrally buoyant plume height, maximum plume induced turbulent diffusivity, and total plume vertically transported water mass flux. These scaling relationships can be used for constructing simplified 1-dimensional models of geochemistry and microbial activity in hydrothermal plumes. Simulation results show that the classical entrainment assumptions, typically invoked to describe hydrothermal plume transport, only apply up to the vertical level of
0.6 times the maximum plume rise height. Below that level, the entrainment coefficient remains relatively constant (
0.15). Above that level, the plume flow consists of a pronounced lateral spreading flow, two branches of inward flow immediately above and below the lateral spreading, and recirculation flanking the plume cap region. Both turbulent kinetic energy and turbulence dissipation rate reach their maximum near the vent; however, turbulent viscosity attains its maximum near the plume top, indicating strong turbulent mixing in that region. The parametric study shows that near vent physical conditions, including chimney height and fluid exit velocity, influence plume mixing from the vent orifice to a distance of
10 times the vent orifice diameter. Thus, physical parameters place a strong kinetic constraint on the chemical reactions occurring in the initial particle-forming zone of hydrothermal plumes. & 2014 Elsevier Ltd. All rights reserved.

Keywords: Hydrothermal plume hydrodynamics Hydrothermal plume turbulence Entrainment Turbulent mixing Computational fluid dynamics Turbulence modeling

1. Introduction Hydrothermal cooling of the ocean ridge system results in the transfer of heat and chemicals from the lithosphere to the ocean (Kadko, 1993; Elderfield and Schultz, 1996; Baker, 2007). The plumes produced by high temperature hydrothermal venting transport significant fractions of this heat and chemicals and interact extensively with seawater (Elderfield and Schultz, 1996). These hydrothermal plumes rise high (100 s of m) above the seafloor (Baker et al., 1995). They are rich in dissolved and particulate inorganic phases and create

n

Corresponding author. Tel.: þ 1 508 289 3641; fax: þ 1 508 457 2194. E-mail address: [email protected] (H. Jiang).

http://dx.doi.org/10.1016/j.dsr.2014.06.006 0967-0637/& 2014 Elsevier Ltd. All rights reserved.

regions of steep redox gradients that can be exploited by deep-sea chemosynthetic organisms. Numerous studies indicate multiple couplings between abiotic and biotic plume processes including chemosynthetic microbial stimulation (Cowen et al., 1986; de Angelis et al., 1993; Lesniewski et al., 2012), dissolved metal complexation (Sander et al., 2007; Bennett et al., 2008; Li et al., 2014), and organic and inorganic particle aggregation (Dymond and Roth, 1988; Toner et al., 2009; Breier et al., 2012). The coupling of these biotic and abiotic processes has implications for the transport of hydrothermal material; specifically, the extent to which hydrothermal material settles to the seafloor or mixes into the ocean (Bennett et al., 2008; Tagliabue et al., 2010; Breier et al., 2012). The mechanisms by which these couplings occur and the composition and magnitudes of the fluxes involved remain areas of active research.

42

H. Jiang, J.A. Breier / Deep-Sea Research I 92 (2014) 41–55

Fluid physics within the rising portion of hydrothermal plumes strongly influences these topics. It places constraints on the kinetics of transport and mixing, and consequently, the occurrence and rate of speciation, particle-forming, and microbially-mediated chemical reactions. Hydrothermal plume mixing and transport has been described by (i) analytical solutions to simplified conservation equations for momentum, mass, and heat, and (ii) numerical solutions to the Navier–Stokes equations for fluid motion. Analytical solutions describing plume mixing and transport are based on the classical fluid dynamics theory for one-dimensional time-averaged plume behavior (Morton et al., 1956; Turner, 1973). The theory is based mainly on two assumptions: First, the shapes of mean vertical velocity and mean buoyancy force are similar at all plume heights. Second, the entrainment velocity (Ue) confined in the horizontal direction at the plume edge, which characterizes the rate of entrainment of ambient fluid into the plume, is proportional to the mean plume vertical velocity (W) at that height. Thus, a constant entrainment coefficient, αe, is defined from Ue ¼ αe  W. This second assumption is the celebrated “entrainment assumption”. Based on these assumptions, simplified conservation equations for momentum, mass and buoyancy are solved analytically for vertical profiles of plume width, vertical velocity and buoyancy. In a stratified fluid, the scaling of maximum plume rise height (Zmax) has been given by Z max ¼ 3:76


 Bexit 1=4 N3

;

ð1Þ

where Bexit is the source buoyancy flux and N is the ambient buoyancy frequency. The scaling coefficient was determined to be 3.76 by analyzing literature data of observed oil fires in atmosphere and of laboratory plume experiments (Briggs, 1969), but its variability was also noted (Chen and Rodi, 1980). The derivation of this scaling equation also assumed linear mixing behavior and small local variation of density throughout the plume. However, these assumptions are violated by the hydrothermal plume in the immediate vicinity of the vent. This is because the temperature of hydrothermal fluids is much higher than that of the ambient seawater and the equations of state of the high-temperature hydrothermal fluids are nonlinear. Thus, the actual source buoyancy flux (i.e. that calculated using high-temperature fluid density at the vent orifice) cannot be used in Eq. (1). For Eq. (1) to be applicable, it is crucial to use the thermal expansion coefficient of the ambient seawater to calculate a so-called asymptotic buoyancy flux and then use it in Eq. (1) (Turner and Campbell, 1987; Woods, 2010). Moreover, the entrainment coefficient is not constant (Fischer et al., 1979; Chen and Rodi, 1980; Wang and Law, 2002). Nevertheless, the theory has been successfully applied to many hydrothermal plume studies and has produced results that compare favorably with observations (e.g. Converse et al., 1984; Lupton et al., 1985; Middleton and Thomson, 1986; Little et al., 1987; Baker et al., 1989; Cann and Strens, 1989; Speer and Rona, 1989; McDougall, 1990; Rudnicki and Elderfield, 1992; Kim et al., 1994; Stranne et al., 2010). The Navier–Stokes equation-based numerical approach allows more detailed prediction than the analytical approach permits. It was used to provide information for the plume cap region and exterior recirculation zone where the entrainment assumption is not valid (Lavelle, 1994). It was used to account for two and threedimensional variations of the plume, such as plume lateral variations (Lavelle, 1994) and spatial variations caused by source spatial distributions (Lavelle, 1995) and by background flow and rotation (Lavelle, 1997; Lavelle et al., 2013). It was also used to study effects of rotation on plume evolution at a relative larger spatial scale and longer time scale (Speer and Marshall, 1995). Time-dependent implementation of the approach also allowed investigating the

initial time development and maturation of a plume caused by episodic hydrothermal discharge (Lavelle and Baker, 1994). Specifically, Lavelle (1997) solved numerically the non-hydrostatic Navier–Stokes equations with the Boussinesq approximation and rotation, together with the incompressible continuity equation and the heat and salt conservation equations. Background flow condition and background profiles of temperature and salinity were properly constructed via solving the governing equations in the steady state without heat forcing and under suitable boundary conditions. Turbulent mixing (both horizontal and vertical) was parameterized as the combination of a constant background term and a term represented by an isotropic mixing coefficient. The isotropic mixing coefficient was explicitly dependent on fluid shear, shear Richardson number, turbulent Prandtl number, turbulence length scale and Smagorinsky constant. The equation of state for seawater (i.e. in the low temperature range; Fofonoff and Millard, 1983) was used for fluid properties. An even-spacing Cartesian mesh was used with a 5 m horizontal spacing and a slightly different vertical spacing. The vent condition was modeled by a given heat source spanning two grid cells at the bottom. Although these numerical models did not explicitly include the individual vent conditions (i.e. exit fluid temperature, exit fluid velocity, and vent geometry) and heavily parameterized turbulent mixing, they generated important insights on the effects of background flow, rotation, and three-dimensionality on hydrothermal plume dynamics and pointed out the importance of turbulent mixing in shaping plume development. Recently, Tao et al. (2013) used a turbulent viscosity (μt) to model hydrothermal plume turbulent mixing and entrainment more sophisticatedly. The evolution equations of both turbulent eddy size (l) and turbulent kinetic energy (k) were solved to determine μt at each grid point and time. Non-uniform Cartesian mesh was used with a 0.59 m finest resolution. The vent condition was explicitly considered by a velocity inlet condition with prescribed high fluid temperature and exit fluid velocity. However, the fluid density was set to be linearly related to temperature. The focuses of all these previous numerical modeling efforts were not on investigating the effects of realistic vent condition and geometry on near vent plume development. Thus, details of model results near the plume stem had to be regarded with caution (Lavelle and Baker, 1994). The focus of this study is on the near vent plume region because this is where some of the most significant changes in plume chemical composition occur. An important goal of this study is the derivation of scaling relationships that can be used for constructing simplified 1-dimensional models of geochemistry and microbial activity in hydrothermal plumes. The approach used is the generation of a numerical simulation that matches a combination of measured vent geometrical, temperature and fluid velocity conditions, observed by remotely operated vehicle, and measured near vent background seawater temperature and stratification, observed by CTD casts. The simulation is then used to provide estimates of those vent conditions and physical driving forces that are not directly measurable. Specifically, we use a computational fluid dynamics (CFD) approach to model a hydrothermal plume from the Eastern Lau Spreading Center: A1 vent in the ABE hydrothermal field (Mottl et al., 2011). We employ the k–ε turbulent model (where k is turbulent kinetic energy and ε is turbulence dissipation rate) to represent turbulent mixing and entrainment. We test two different k–ε turbulence closure models for their suitability for hydrothermal plume modeling. We construct a non-uniform triangular mesh to explicitly represent vent geometry. We apply a set of nonlinear equations of state for high-temperature hydrothermal fluids and seawater to represent vent fluid properties as realistically as possible in the near vent region (Sun et al., 2008). We model a baseline plume case using measurements and

H. Jiang, J.A. Breier / Deep-Sea Research I 92 (2014) 41–55

z

estimates of near vent parameters. We conduct a series of parametric sensitivity studies to (i) generate numerical data for calculating the entrainment coefficient (αe) as a function of the vertical distance above the vent and for deriving the scaling of turbulent diffusivity within the rising plume, (ii) compare the performances of the two turbulence closure models, and (iii) determine the relative importance of vent temperature, vent exit velocity, background stratification, and chimney height on plume mixing and transport. We validate the model with comparisons to water column observations of plume structure from the ABE-A1 site. We investigate the effect of near vent physical conditions, including chimney height and fluid exit velocity, on plume entrainment and mixing. We derive scaling relationships for simplified 1-dimensional biogeochemical models concerning plume impact on the ocean.

2.1. Computational parameters Our plume flow model is axisymmetric with the axis at the vent centerline (Fig. 1). The horizontal span (i.e. in the radial direction) of the computational domain is 200 m. The depth (i.e. in the axial direction) of the domain is 400 m, which is about twice the plume rise height achieved by our baseline modeling case. This is similar to the choices of domain size made by previous plume flow modeling studies (Lavelle, 1997; Tao et al., 2013). Suzuki and Koyaguchi (2010) have conducted three-dimensional numerical simulations of the flow of terrestrial volcanic eruption clouds. Their simulations have suggested that the mesh size in the near vent region should be sufficiently smaller than D/10 (where D is vent diameter) in order to correctly reproduce turbulent mixing near the vent. Our choice of the finest mesh size is based on their suggestion; in specific, the vent radius (¼ 0.07 m) is divided evenly into 10 segments. Thus, the smallest mesh size is 0.007 m. The edge at the top of the domain is discretized evenly into 50 segments, leading to a mesh size of 4 m at the domain top. The whole domain is discretized by
56,000 triangular cells that are

Pressure outlet Constant Ttop

Axisymmetric axis

400 m

2. Methods

ZOOM

Wexit, Texit

Our hydrothermal plume flow modeling was performed with the finite-volume CFD software package ANSYS FLUENT (version 13.0.0). The governing equations are the density weighted unsteady Reynolds-averaged Navier–Stokes (URANS) equations and energy equation, which are numerically solved for the mean velocity and temperature fields. The effects of turbulence are represented by the so-called Reynolds stresses that appear as unknowns in the URANS equations. The URANS equations are closed by relating the Reynolds stresses to the mean velocity gradients via the turbulent viscosity hypothesis. Analogically, the Reynolds-averaged energy equation is closed by relating the temperature fluxes (as unknowns) to the mean temperature gradients via the gradient-diffusion hypothesis. In this study, we employ a two-equation k–ε turbulence model for the closure, where two transport equations are solved for turbulent kinetic energy (k) and the rate of dissipation of turbulent kinetic energy (ε), respectively. We apply and compare two k–ε turbulence models, i.e. the realizable k–ε model and the renormalization group (RNG) k–ε model. Both models include the effects of buoyancy on the production and dissipation of turbulent kinetic energy. These turbulence modeling equations are rather generic; we therefore present them in a separate document (equations_numerics.pdf) supplied as electronic supplementary material. We also include details of numerical methods in the same document. In the following, we present model setup and calculation methods that are directly related to plume flow modeling.

43

No-slip 200 m Constant Tbottom

r

H=2 m

Exit diameter D=0.14 m

Fig. 1. Grid and boundary conditions of an axisymmetric computational fluid dynamics model for simulating the ABE-A1 vent hydrothermal plume. Upper panel: the whole computational domain; Lower panel: the near vent region.

stretched at a constant rate of 1.02 from the vent orifice to the domain boundaries (Fig. 1 upper panel). The resulting mesh sizes are less than 1 m within the bottom 100 m of the domain. The geometry of the vent chimney is represented by the triangular mesh (Fig. 1 lower panel). At the vent orifice, the velocity-inlet boundary condition is prescribed with constant exit fluid velocity (Wexit) and temperature (Texit ¼309 1C). At the bottom of the domain, i.e. the seafloor, a no-slip boundary condition is prescribed with a constant temperature (Tbottom ¼ 2.2 1C). At the top of the domain, a pressure-outlet boundary condition is prescribed with a constant temperature (Ttop, depending on the background stratification different cases may have different Ttop'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 initial condition for the flow field is zero flow velocity everywhere except at the vent orifice where the vertical flow velocity is Wexit. The initial condition for the temperature field is a vertically linear varying temperature field everywhere from Tbottom at the seafloor to Ttop at the domain top, except at the vent orifice where the temperature is Texit( ¼309 1C). The vertically linear varying temperature field determines the background buoyancy frequency (N), which is calculated according to: N¼

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi g ρðT bottom Þ  ρðT top Þ ; H top  H ρðT bottom Þ

where For where use a

ð2Þ

all the symbols are described in Table 1. each simulation case, the total integration time is 4tn, tn is the buoyancy time scale, defined as t n ¼ 2π =N. We second order implicit time integration scheme in our

44

H. Jiang, J.A. Breier / Deep-Sea Research I 92 (2014) 41–55

Table 1 Symbols used and default parameter values. Symbol Explanation

Units

Default value

g T ρ β Cp μ κ D H Htop Ttop Tbottom Texit N

Gravitational acceleration Temperature Mass density of seawater Thermal expansion coefficient of seawater Specific heat of seawater Dynamic viscosity of seawater Thermal conductivity of seawater Effective diameter of vent orifice Vent chimney height above seafloor Height of the computational domain Temperature at the model top (i.e. Htop above seafloor) Temperature at seafloor Temperature of vent exit fluid qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ρðT bottom Þ  ρðT top Þ g Background buoyancy frequency, calculated as N ¼ ρðT bottom H top  H Þ

m s2 1C kg m  3 K1 J kg  1 K  1 kg m  1 s  1 W m1 K1 m m m 1C 1C 1C s1

9.81 – – – – – – 0.14 – 400 – 2.2 309 –

Wexit Iexit

Upward velocity of vent exit fluid Turbulence intensity (defined as the ratio of the root-mean-square of the velocity fluctuations to the mean flow velocity) at vent orifice

m s1 %

– –

Vexit

Volume flux issued from vent orifice, calculated as V exit ¼ W exit πD4

m3 s  1

–

kg s  1

–

kg s  1 W m4 s  3

– – –

m m

– –

– m s1 – J kg  1 W kg  1 m2 s  1 –

– – – – – – –

Mexit Mtotal Qexit Bexit Zmax Zneutral tn w C k ε νt αe

Fluid mass flux issued from vent orifice, calculated as M exit ¼ Total plume vertically transported water mass flux Heat flux issued from vent orifice, calculated from CFD

2

2 ρðT exit ÞW exit πD4

Þ  ρðT exit Þ Buoyancy flux issued from vent orifice, calculated as Bexit ¼ g ρðT bottom V exit ρðT bottom Þ

Maximum plume rise height where the plume centerline vertical velocity w becomes 0 for the first time Neutrally buoyant height of the plume where the plume centerline density becomes equal to the density of the background seawater for the first time Buoyancy time scale, defined as t n ¼ 2π=N Vertical velocity component Scalar concentration normalized by the scalar concentration of the vent exit fluid (i.e. C¼ 1 at the vent orifice) Turbulent kinetic energy Turbulence dissipation rate (Kinematic) turbulent viscosity Entrainment coefficient

simulations. An implicit scheme is unconditionally stable with respect to time step size. However, we choose a small time step of 0.05 s to ensure that the solution is converged within 20 iterations at each time step. Because the total integration time is relatively short, we neglect Earth’s rotation in our plume flow modeling. To illustrate the transport and dilution of vent fluid derived chemical constituents we also include a chemically non-reactive scalar quantity, C, in our model. This is a generic scalar quantity but the distribution of this scalar within the model should be representative of the relative distribution of conservative vent fluid chemical constituents within the plume. Conceptually, helium-3 is one constituent that should follow this distribution pattern and thus we consider this scalar to represent helium-3. The scalar distribution is numerically solved from a steady-state scalar transport equation. The steady-state plume flow field (i.e. the simulated flow field at 4tn) is used in the advection term of the scalar transport equation. For the diffusion term of the scalar transport equation, the effective scalar diffusivity, deff, is determined as deff ¼d þmt/Sct, where d is the molecular diffusivity, mt is the turbulent viscosity, and Sct is the turbulent Schmidt number (¼ 0.7). The boundary conditions for C are C ¼1 at the vent orifice, C ¼0 at the top boundary, and freely adjustable at the lateral boundary.

2.2. Method for calculating the entrainment coefficient from numerical simulation data Using our numerical data, for each case we directly determine the local (i.e. z-dependent) entrainment coefficient (αe) from the following equation (e.g. Pham et al., 2006; Plourde et al., 2008;

Suzuki and Koyaguchi, 2010): dM ae ¼ dz ; Ω

ð3Þ R1

where MðzÞ  0 2πρwrdr is the vertical mass flux inside the plume [where z is the axial (i.e. vertical) coordinate, r is the radial (i.e. horizontal) coordinate, ρ(z, r) is time-averaged density at each grid point, and w(z, r) is time-averaged vertical velocity at each qRffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 2 rdr [where ρ (z) is the grid point], and ΩðzÞ  2πρ0 ðzÞ w 0 0 background density]. Our numerical simulations are timedependent. The simulated data are output at every 60 s over the total integration time of 4tn. ρ(z, r) and w(z, r) are obtained by time-averaging those time series data that are output during the last tn for each case. 2.3. Parametric sensitivity studies Our parametric sensitivity studies vary parameters at two levels. At the highest level, two different turbulence closure models, the realizable k–ε model (rke) and the RNG k–ε model (rngke), are employed in two separate groups of simulations (Tables 2 and 3). These datasets are used to evaluate the performance of these two turbulence models for predicting hydrothermal plume physical transport. Within the two turbulence model groups, different values for background buoyancy frequency (N) are achieved by setting different temperatures for Ttop (i.e. rke1-6 and rke7_0 in Table 2, and rngke1-6 and rngke7_0 in Table 3). With all other parameters the same, simulation runs using different values of N generate datasets with different maximum plume rise height (Zmax). These

H. Jiang, J.A. Breier / Deep-Sea Research I 92 (2014) 41–55

45

Table 2 Summary of the simulation runs using the realizable k–ε model (rke). Simulation run

H (m)

Tbottom (1C)

Ttop (1C)

Wexit (m s  1)

Iexit (%)

Vexit (m3 s  1)

N (s  1)

rke1 rke2 rke3 rke4 rke5 rke6 rke7_0n rke7_1 rke7_2 rke7_3

2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 0.0

2.2 2.2 2.2 2.2 2.2 2.2 2.2 2.2 2.2 2.2

6.2 5.2 4.2 3.7 3.2 2.7 2.3663 2.3663 2.3663 2.3663

0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.1 0.2 0.2

0 0 0 0 0 0 0 0 10 0

0.00307876 0.00307876 0.00307876 0.00307876 0.00307876 0.00307876 0.00307876 0.00153938 0.00307876 0.00307876

0.002982 2662,117.4 0.00881905 0.002494 2662,117.4 0.00881905 0.001884 2662,117.4 0.00881905 0.001574 2662,117.4 0.00881905 0.001258 2662,117.4 0.00881905 0.000874 2662,117.4 0.00881905 0.000473 2662,117.4 0.00881904 0.000473 1331,627.4 0.00440952 0.000473 2667,914.8 0.00881904 0.000472 2663,528.9 0.00881904

n

Qexit (W)

Bexit (m4 s  3)

Zmax (m)

Zneutral (m)

Mexit (kg s  1)

Mtotal (kg s  1)

58.469 65.013 75.558 84.165 97.900 127.702 201.255 168.353 203.724 202.979

42.670 47.714 55.477 61.849 72.311 94.586 151.397 125.536 152.647 152.424

2.25908 2.25908 2.25908 2.25908 2.25908 2.25908 2.25908 1.12954 – –

4,314.28 5,159.98 6,564.99 7,902.60 10,237.51 15,745.73 30,306.33 18,627.54 – –

Baseline modeling case, which has been designed to approximately reproduce the observed plume flow issued from the ABE-A1 vent (see Section 2.5).

Table 3 Summary of the simulation runs using the RNG k–ε model (rngke). Simulation run

H (m)

Tbottom (1C)

Ttop (1C)

Wexit (m s  1)

Iexit (%)

Vexit (m3 s  1)

N (s  1)

Qexit (W)

Bexit (m4 s  3)

Zmax (m)

Zneutral (m)

rngke1 rngke2 rngke3 rngke4 rngke5 rngke6 rngke7_0 rngke7_1 rngke7_2 rngke7_3

2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 0.0

2.2 2.2 2.2 2.2 2.2 2.2 2.2 2.2 2.2 2.2

6.2 5.2 4.2 3.7 3.2 2.7 2.3663 2.3663 2.3663 2.3663

0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.1 0.2 0.2

0 0 0 0 0 0 0 0 10 0

0.00307876 0.00307876 0.00307876 0.00307876 0.00307876 0.00307876 0.00307876 0.00153938 0.00307876 0.00307876

0.002982 0.002494 0.001884 0.001574 0.001258 0.000874 0.000473 0.000473 0.000473 0.000472

2662,112.5 2662,108.7 2662,111.5 2662,110.7 2662,111.4 2662,111.6 2662,112.8 1373,623.5 2640,131.4 2663,451.3

0.00881905 0.00881905 0.00881905 0.00881905 0.00881905 0.00881905 0.00881904 0.00440952 0.00881904 0.00881904

55.808 62.526 73.268 81.930 96.120 126.286 195.196 164.668 194.888 196.325

39.407 44.443 52.199 58.577 68.996 91.238 142.552 118.935 142.231 143.593

datasets are used to calculate the entrainment coefficient (αe) as a function of the non-dimensional vertical distance above the vent orifice (i.e. the vertical distance normalized by Zmax). Also within the two turbulence model groups (i.e. Tables 2 and 3), three near vent physical parameters are varied. The three parameters are: (i) vent chimney height (rke7_0 and rke7_3 in Table 2; rngke7_0 and rngke7_3 in Table 3), (ii) vent exit fluid velocity (rke7_0 and rke7_1 in Table 2; rngke7_0 and rngke7_1 in Table 3) and (iii) vent exit fluid turbulent condition (rke7_0 and rke7_2 in Table 2; rngke7_0 and rngke7_2 in Table 3). We use these datasets to investigate the effects of near vent conditions on the entrainment coefficient. We use the datasets that vary actual vent buoyancy flux (Bexit) and the background buoyancy frequency (N) to determine the scaling of Zmax about Bexit and N: Z max ¼ A1 ðBexit N  3 Þ1=4 :

ð4Þ

As discussed previously, the scaling coefficient A1 in this case will not be equal to 3.76 because here Bexit is not the so-called asymptotic buoyancy flux but the actual vent buoyancy flux: Bexit ¼ g

ρðT bottom Þ  ρðT exit Þ V exit ; ρðT bottom Þ

ð5Þ

where ρ(Tbottom) is the density at Tbottom, ρ(Texit) is the density at Texit, and Vexit is volume flux at the vent orifice (Table 1). We also use the datasets to determine the scaling of total plume vertically transported water mass flux (Mtotal) with Bexit and N: M total ¼ A2  ρðT bottom Þ  ðBexit 3 N  5 Þ1=4

ð6Þ

where A2 is the scaling coefficient. Finally, we use those datasets that vary Bexit and N to determine the scaling of turbulent viscosity along the vent centerline about Bexit and N.

2.4. Hydrothermal fluid, seawater, and vent properties The density (ρ) and specific heat (Cp) of hydrothermal fluid and seawater are calculated from a set of nonlinear equations of state of a thermal saline fluid for a temperature range of 0–374 1C, pressure range of 0.1–100 MPa and absolute salinity range of 0–40 g kg  1 (Sun et al., 2008; Fig. 2A and B). To our knowledge, there are no similar equations available for calculating thermal conductivity and viscosity of the thermal saline fluid in the similar temperature, pressure and salinity ranges. We therefore calculate thermal conductivity and viscosity using equations available for pure water according to National Institute of Standards and Technology (NIST) Reference Fluid Thermodynamic and Transport Properties Database (REFPROP), NIST Standard Reference Database 23, Version 9.0 (Fig. 2C and D). Because the plume flow is turbulent, turbulent viscosity and thermal diffusivity will dominate over fluid viscosity and thermal conductivity, respectively. Thus, the choice of pure water viscosity and thermal conductivity should not affect our simulation results. Within the rise height of the plume (
200 m), pressure change has negligible effects on the fluid properties compared to those caused by temperature change (Fig. 2). Though salinity may vary near the vent, in situ CTD measurements indicate that for much of the rising plume salinity only varies within a very small range around 34.65. Thus, we calculate water properties only as functions of temperature at fixed salinity (S ¼34.65) and pressure (¼ 21.56 MPa, i.e. the vent water depth of 2140 m). 2.5. Baseline modeling case The goal of our baseline modeling case is to provide a best estimate of steady-state plume flow, in the absence of bottom currents, from a vent (i.e. A1, 201450 47.5″S 1761110 28.8″W, depth 2140 m) in the northern region of the ABE hydrothermal field along the Eastern Lau Spreading Center (Ferrini et al., 2008; Mottl et al., 2011). Note, there are many vents in this field. There are at

46

H. Jiang, J.A. Breier / Deep-Sea Research I 92 (2014) 41–55

background

Seawater (S=34.65) at 21.56 MPa (water depth ~2140 m) Seawater (S=34.65) at 20.15 MPa (water depth ~2000 m)

Pure water at 21.56 MPa (water depth ~2140 m) Pure water at 20.15 MPa (water depth ~2000 m) background

Fig. 2. Seawater properties used in simulations, (A) mass density, ρ, (B) specific heat, Cp, (C) thermal conductivity, κ, and (D) dynamic viscosity, m, as functions of temperature, T.

least 8 areas of localized high and low temperature hydrothermal activity within the ABE field, and many of these areas contain multiple high temperature vent edifices (Ferrini et al. 2008). This particular region of the field and vent were selected for modeling because we had supporting measurements, collected by remotely operated vehicle (ROV) and CTD profiles, from (i) the rising plume emanating from A1 vent and (ii) the neutrally buoyant plume overlying this region of the ABE hydrothermal field. A1 vent is a typical high temperature vent within the ABE field, but it does not necessarily represent the average parameters of vents within the ABE field. 2.6. Field measurements The effective vent orifice diameter was estimated based on a video imagery collected by ROV Jason II (cruise TN235, Dive J2-424). Vertical plume and background seawater structure with respect to temperature, density and optical backscatter was identified by multiple nearby CTD profiles, on a water rosette equipped with a Seapoint turbidity sensor, during R/V Thomas G. Thompson cruises TN235 and TN236. Rising plume temperature

Fig. 3. CTD profile data from field experiments conducted around the ABE-A1 plume. (A) Vertical profiles of potential temperature. (B) Vertical profiles of potential density. (C) Three plume CTD profiles with temperature anomalies compared to the background density profile. Note that these CTD profiles passed through the non-buoyant plume but did not intersect the buoyant plume.

measurements were made with the temperature probe of the In situ Electrochemical Analyzer system (Luther et al., 2008). The maximum measured temperatures from 30 to 60 s intervals at each plume elevation were reported.

3. Results Within the local environment of the ABE hydrothermal field, the near vent background seawater temperature was
2.2 1C based on available in situ CTD measurements. The CTD measured temperature profiles suggest a neutrally buoyant plume height (Zneutral) of
150 m (Fig. 3A). The CTD measurements also reveal a density anomaly with respect to background seawater up to an elevation of 200 m (Fig. 3B and C), suggesting a maximum plume rise height (Zmax) of
200 m. From the CTD data, the background buoyancy frequency (N) is estimated to be
0.0005 s  1. The exit fluid vertical velocity (Wexit) could not be measured directly. We estimated Wexit by iteratively varying this parameter to numerically generate a

H. Jiang, J.A. Breier / Deep-Sea Research I 92 (2014) 41–55

w (m s-1)

2.0 1.5 (z-H)/Zmax

plume of Zmax
200 m. By this method our best estimate of Wexit is 0.2 m s  1. This is generally consistent with typical hydrothermal vent discharge velocities, which range between 0.1 and 2 m s  1 (Ramondenc et al., 2006). The ABE-A1 vent chimney structure consisted of at least four orifices in close proximity; each orifice was
0.07 m in diameter. The total area of the four orifices was
0.0154 m2, equivalent to an effective orifice diameter of
0.14 m. Thus, we modeled the ABE-A1 vent as having a single orifice of 0.14 m in diameter. The maximum exit fluid temperature measured by Mottl et al. (2011) for this vent was 309 1C.

47

1.0 0.5 0.0

3.1. Baseline modeling case: ABE-A1 vent

k (J kg-1)

2.0

3.1.1. Temporal evolution of model plume Temporal evolution of flow velocity and temperature along the vent centerline (Fig. 4A) shows that the initial plume overshoots its steady state maximum rise height (Zmax) before one tn (i.e. the buoyancy time scale), and then oscillates around Zmax for several tn’s before finally reaching steady state. Temporal evolution of turbulent kinetic energy (k, Fig. 4B) and turbulence dissipation rate (ε, Fig. 4C) shows similar trends. Thus, a total integration time of 4tn should be sufficient to achieve the steady state of the model plume. To visualize the temporal evolution of the model plume, we have provided animations of flow velocity (rke7_0_velVector. mp4), temperature (rke7_0_temp.mp4), turbulent kinetic energy (rke7_0_tke.mp4), turbulence dissipation rate (rke7_0_epsilon. mp4), and turbulent viscosity (rke7_0_turbVisco.mp4) as electronic Supplementary material. 3.1.2. Steady-state model plume Modeled plume vertical velocity reaches its maximum (
1.27 m s  1) just
0.30 m above the vent orifice (Fig. 5A). This maximum velocity is more than 6 times the plume exit velocity (Wexit ¼0.2 m s  1). Such a rapid acceleration over a very short vertical distance is caused by the extremely high buoyancy of the exit fluid. Above this short vertical distance, the plume vertical velocity starts to decrease due to seawater entrainment and subsequent buoyancy loss. The height where the plume centerline vertical velocity initially reaches zero is the maximum plume rise height (i.e. Zmax). Modeled plume centerline temperature decreases from 309 1C at the vent orifice to
167 1C at 0.5 m above the vent orifice, to
4.45 1C at 10 m above the vent orifice, and to
2.48 1C at 40 m above the vent orifice (Fig. 5B). The height where the plume centerline temperature initially reaches background seawater temperature is the neutrally buoyant plume height (i.e. Zneutral). In situ temperature measurements made within the ABE-A1 plume, specifically at 10 and 40 m above the vent orifice, are consistent with the model predictions (Fig. 5B). Because the horizontal temperature gradient near the vent orifice is very high, temperature measurements made at locations slightly off the vent centerline in this region will be much less than the centerline

(z-H)/Zmax

1.5 1.0 0.5 0.0 ε (W kg-1)

2.0 (z-H)/Zmax

Simulation run rke7_0 (Table 2, realizable k–ε model) is our best estimate of physical transport and mixing for a contiguous hydrothermal plume representing both the ABE-A1 rising plume and the CTD observations of the overlying neutrally buoyant plume in this region on the ABE field. This is our baseline modeling case. It is worth noting that while it is based on measurements and estimates of actual vent conditions it is still a simplified representation of the system; the potential limitations of these simplifications are discussed in Section 4. This section describes the temporal evolution of this model, the steady state conditions, turbulence within the plume, and plume structure with respect to temperature, flow, and distribution of a model, inert, vent-fluid enriched, chemical constituent representing helium-3.

1.5 1.0 0.5 0.0 0

1

2 t/t*

3

4

Fig. 4. Results from the baseline modeling case (rke7_0, Table 2). Time evolutions of vertical profiles along the vent centerline of a starting plume, for (A) vertical velocity, w (filled color contours) overlapped with temperature (in 1C, black line contours), (B) turbulent kinetic energy, k, and (C) turbulence dissipation rate, ε. z is vertical coordinate, H is the vent orifice height above the seafloor, Zmax is maximum height of rise of the plume, t is time, and t* is the buoyancy time scale defined as t* ¼2π/N where N is the ambient buoyancy frequency.

temperature at the same height. The in situ temperature measurement of 7 1C at 0.5 m above the ABE-A1 vent does indeed fall well below the predicted centerline temperature. In addition, actual plume temperatures may well have varied temporally due to timedependent background currents not considered by our simulations. A more realistic comparison of modeled and measured plume temperatures at the 0.5 m plume elevations involves some horizontal averaging of the predicted temperatures across the plume profile. Averaged horizontally across the 20 cm full-width of the plume, the modeled temperature at 0.5 m above the vent is
80 1C. This is still approximately 10 times greater than the measured temperature. This suggests that either (i) the temperature measurement at 0.5 m was still further removed from the vent centerline and/or (ii) that some of our simplifying model assumptions and estimated parameters are not fully representative of the conditions at the time of these measurements; these points are discussed in Section 4. However, modeled centerline temperatures should represent upper limits for measured temperatures at these plume heights for a single contiguous plume that reaches a neutrally buoyant plume elevation of 150 m. The modeled scalar, representing helium-3, decreases from 1 at the vent orifice to
0.8354 at 0.5 m above the vent orifice, to
0.0820 at 10 m above the vent orifice, and to
0.0225 at 40 m above the vent orifice (Fig. 5C).

48

H. Jiang, J.A. Breier / Deep-Sea Research I 92 (2014) 41–55

CFD simulated, along axis

CFD simulated, along axis background in CFD

CFD simulated, along axis

Zmax

Zmax Zneutral Temperature measurements: 0.5 m above vent 10 m above vent 40 m above vent

Fig. 5. Results from the baseline modeling case (rke7_0, Table 2). (A) Vertical velocity, w, profile along the vent centerline at the time instant 4t* after the start of the plume. (B) Temperature, T, profile along the vent centerline at the time instant 4t* after the start of the plume. Also plotted are background temperature profile and temperature measurements made around the ABE-A1 plume at three different heights above the vent orifice. (C) Steady state scalar concentration, C, profile along the vent centerline, where C¼ 1 at the vent orifice and C¼0 at the top boundary.

300 k (J kg-1)

ε (W kg-1)

ν t (m2s-1)

250 Zmax

200

backg

round

150

98

100 50 0 97

Fig. 6. Results from the baseline modeling case (rke7_0, Table 2). Contours of (A) turbulent kinetic energy, k, (B) turbulence dissipation rate, ε, and (C) turbulent viscosity, νt, plotted for the time instant 4t* after the start of the plume. Also plotted in (D) are four beam transmission vertical profiles from field experiments conducted around the ABEA1 plume.

3.1.3. Modeled hydrothermal plume turbulence Hydrothermal plume flow is highly turbulent. Contour plots of turbulent kinetic energy (k, Fig. 6A), turbulence dissipation rate (ε, Fig. 6B) and turbulent viscosity (νt ¼ μt/ρ, Fig. 6C) all display a normally expected mushroom-like spatial distribution, with elevated values confined inside a cone-shaped region emanating from the vent orifice until
120 m and a cap-shaped region between
120 and 200 m (i.e. the maximum plume rise height). Along the vent centerline, turbulent kinetic energy (k) reaches its maximum, i.e.
0.045 J kg  1, at
0.65 m above the vent orifice and decreases gradually to
8.3  10  5 J kg  1 at
200 m above the vent (Fig. 6A); turbulence dissipation rate (ε) reaches its maximum, i.e.
0.13 W kg  1, at
0.57 m above the vent orifice and decreases rather rapidly to
1.9  10  8 W kg  1 at
200 m above the vent (Fig. 6B); turbulent viscosity (νt) is above 1.0  10  4 m2 s  1 from
0.21 to
222 m above the vent and reaches its maximum, i.e.
6.85  10  2 m2 s  1, at
186 m above the vent orifice (Fig. 6C). Values of oceanic turbulence dissipation rate range considerably from
10  10 W kg  1 in the abyssal ocean to
10  1 W kg  1 in surf zones, straits and regions with very rapid tidal currents (Thorpe, 2005). Thus, within the entire rise height of the plume, the values of plume turbulence dissipation rate are at least two orders of magnitude higher than the typical values in the abyssal ocean. Along the vent centerline from
0.34 to
0.97 m above the vent orifice, the values of plume turbulence dissipation rate are all above 0.05 W kg  1, which are comparable to the values of most energetic regions in the ocean. Thurnherr and St. Laurent (2012) reported a microstructure profiler and CTD measured patch of strong mixing with peak dissipations exceeding 1.5  10  8 W kg  1 associated

with a buoyant hydrothermal plume. Our simulated dissipation at the plume maximum rise height is comparable in magnitude to their observed dissipation. Within the plume flow, turbulent viscosity is not constant but varies greatly in space (Fig. 6C). Note that the values reported here are only for our baseline modeling case and will differ depending on vent conditions, which vary considerably. Among the three plume-intersecting vertical CTD profiles, two (Fig. 6D, pink and red lines) have minimum beam transmissions at elevations coincident with our neutrally buoyant plume height model prediction. This is consistent with an enriched hydrothermal mineral particle laden plume. One CTD profile (Fig. 6D, blue line) shows a minimum beam transmission at a slightly higher elevation. This may be due to the unsteadiness of the plume exit flux, which was not directly considered in our baseline modeling case.

3.1.4. Modeled steady-state plume structure Modeled plume flow shows a persistent entrainment region below the level of
0.6Zmax. This region is characterized by nearly horizontal inward flow at the edge of the plume (Fig. 7A). Related to this, the local entrainment coefficient, αe, has been calculated to be roughly constant (
0.15) for vertical levels below
0.6Zmax (Fig. 7B). The local entrainment coefficient, αe, decreases from
0.15 at
0.6Zmax to 0 at
0.72Zmax. Above 0.72Zmax, αe becomes negative because the flow is dominated by lateral spreading centered around
0.8Zmax. However, above and below this strong lateral spreading flow there are two branches of weak inward flow

H. Jiang, J.A. Breier / Deep-Sea Research I 92 (2014) 41–55

49

Fig. 7. Results from the baseline modeling case (rke7_0, Table 2). (A) Plume flow velocity vector field (colored vectors with uniform vector length; only vectors with magnitude Z 0.001 m s  1 are shown) overlapping with contours of temperature field (red line contours) for the time instant 4t* after the start of the plume. Also overlapped are contours of steady-state scalar concentration (black line contour levels: 3.98  10  3, 4.47  10  3, 5.01  10  3, 5.62  10  3, 6.31  10  3, 7.08  10  3, 7.94  10  3, 8.91  10  3, and 1.0  10  2, relative to 1 at the vent orifice and 0 at the top boundary). The region lateral to the plume core and below Zmax, including the region of the strong lateral spreading flow, has a relatively homogenous scalar distribution, with concentration ranging from 6.31  10  3 to 7.08  10  3. (B) Vertical profile of the local entrainment coefficient, αe.

(Fig. 7A). Along the boundaries between the lateral spreading flow and the two branches of inward flow, shear instabilities are likely as evidenced by flow curvature and eddy formation in these areas. Between
0.6Zmax and Zmax and lateral to the plume core region, there is strong recirculation formed by (i) the upward flow of the plume core, (ii) the lateral spreading flow originating from the plume cap region and (iii) the entrained inward flow at
0.6Zmax. The plume flow pattern above
0.6Zmax deviates greatly from the self-similar entrained flow pattern below
0.6Zmax. Thus the entrainment assumption is not valid above
0.6Zmax and values of αe calculated here for elevations above
0.6Zmax using the simulation data are not consistent with their original physical meaning.

rke

rngke

3.2. Parametric sensitivity The turbulence model used in the plume simulations influences the value of αe, and thus the degree of mixing within the simulations. The RNG k–ε model (rngke, Fig. 8B) generates larger values of αe than does the realizable k–ε turbulence model (rke, Fig. 8A). As a result, the maximum plume rise heights achieved by the RNG k–ε model are smaller than those achieved by the realizable k–ε model (Tables 2 and 3). For the realizable k–ε model (Fig. 8A), αe increases quickly from
0.1 at the bottom level to
0.15 at
0.03Zmax, ranges between 0.14 and 0.16 from
0.03Zmax to
0.6Zmax, and then decreases from
0.14 at
0.6Zmax to 0 at
0.72Zmax. For the RNG k–ε model (Fig. 8B), αe increases quickly from
0.1 at the bottom level to
0.19 at
0.03Zmax, decreases gradually from
0.19 at
0.03Zmax to
0.15 at
0.6Zmax, and then decreases from
0.15 at
0.6Zmax to 0 at
0.72Zmax. Vent exit velocity is inversely related to mixing within the near vent region where elevation is o 10 times the vent orifice diameter. Specifically, the models (both rke and rngke variants) with higher exit flow velocity (e.g. Wexit ¼0.2 m s  1, rke7_0 in Table 2) result in smaller local entrainment coefficients (αe) than models with a lower exit flow velocity (e.g. Wexit ¼0.1 m s  1, rke7_1 in Table 2; Fig. 9A). To compare the relative contributions from buoyancy and momentum for these

Fig. 8. Vertical profile of entrainment coefficient, αe. (A) Determined from eight simulation runs (i.e. rke1-6, rke7_0 and rke7_1) listed in Table 2 using the realizable k–ε model (rke). The triangle symbols are the vertical profiles of αe calculated from these eight simulation runs, and the solid line is the averaged vertical profile of αe. The standard deviation for αe varies vertically with maximum 0.020. (B) Determined from eight simulation runs (i.e. rngke1-6, rngke7_0 and rngke7_1) listed in Table 3 using the RNG k–ε model (rngke). The triangle symbols are the vertical profiles of αe calculated from these eight simulation runs, and the solid line is the averaged vertical profile of αe. The standard deviation for αe varies vertically with maximum 0.025.

two cases, we consider the flux Richardson number, Rf, defined as (see McDuff, 1995): 1=2

Rf 

V exit Bexit ðV exit W exit Þ

5=4

¼

1=2 Bexit ;; 1=4 5=4 V exit W exit

ð7Þ

where Bexit is buoyancy flux, Vexit is volume flux, and Wexit is upward velocity of vent exit fluid. Using values listed in Table 2, we have calculated Rf
3 for rke7_0 and Rf
6 for rke7_1. Since Rf 41 for both cases, the buoyancy contribution dominates over the momentum contribution and thus these fluid motions are best

50

H. Jiang, J.A. Breier / Deep-Sea Research I 92 (2014) 41–55

rke7_0 (Wexit=0.2 m s-1) rke7_1 (Wexit=0.1 m s-1)

rke7_0 (H=2 m)

rke7_0 (Iexit=0)

rke7_3 (H=0)

rke7_2 (Iexit=10 %)

Fig. 9. Effects of different near vent physical conditions (Table 2) on entrainment coefficient, αe, in the near vent region. (A) Different exit velocity of vent fluid: rke7_0 (Wexit ¼ 0.2 m s  1) vs. rke7_1 (Wexit ¼ 0.1 m s  1). (B) Different vent chimney height above seafloor: rke7_0 (H¼ 2 m) vs. rke7_3 (H¼ 0). (C) Different turbulence intensity at vent orifice: rke7_0 (Iexit ¼ 0) vs. rke7_2 (Iexit ¼ 10%). Note that in the plots the vertical distance above the vent orifice, z  H, is normalized by vent orifice diameter, D.

described as buoyant plumes rather than fluid jets. However, rke7_1 has a larger Rf and therefore buoyant flow is even more pronounced with greater mixing in the near vent region. Chimney height is positively related to plume mixing within the near vent region where elevation is o10 times the vent orifice diameter (Fig. 9B). Specifically, models (again both rke and rngke variants) with a taller chimney (e.g. H ¼2 m, rke7_0 in Table 2) have larger αe than models with no chimneys (e.g. H ¼0 m, rke7_3 in Table 2). Turbulence intensity of the vent fluid has no effect on αe (Fig. 9C).

Zmax: Y=1.999X+8.320 2 (R =0.99889)

Zmax: Y=1.958X+7.558 2 (R =0.99954)

Zneutral: Y=1.517X+4.414 2 (R =0.99839)

Zneutral: Y=1.443X+3.796 2 (R =0.99936)

3.3. Maximum plume rise height scaling Using Eqs. (2), (4) and (5), we determine scaling coefficients (A1) for relationships between physical parameters associated with venting and (i) maximum plume rise height (Zmax) and (ii) neutrally buoyant plume height (Zneutral), respectively. In these relationships, the physical parameters associated with venting are vent temperature, vent exit velocity and bottom seawater temperature. These parameters are represented collectively by the buoyancy flux (Bexit) and the background buoyancy frequency (N) terms in the scaling relationship (Eq. (4)). For the realizable k–ε turbulence model (rke), the resulting relationships are Zmax ¼1.999 (Bexit/N3)1/4 and Zneutral ¼ 1.517 (Bexit/N3)1/4 (Fig. 10A). For the RNG k–ε model (rngke), the resulting relationships are Zmax ¼1.958 (Bexit/N3)1/4 and Zneutral ¼ 1.443 (Bexit/N3)1/4 (Fig. 10B). The Zneutral/ Zmax ratios for the rke and rngke models are 0.759 and 0.737, respectively, which are close to the value of 0.761 derived from analytical solutions of plume transport based on classical fluid dynamics theory (Pages 197–198 in Turner, 1973). From a practical standpoint, the derived scaling relationship for the neutrally buoyant plume height (Zneutral) can be used to estimate vent conditions from observation of the neutrally buoyant layer. Alternatively, Turner and Campbell (1987) based similar scaling relationships on the so-called asymptotic buoyancy flux (Basymp) in place of Bexit, where Basymp ¼gβbottom(Texit  Tbottom)Qexit with βbottom the thermal expansion coefficient of bottom seawater (¼ 6.203  10  5 K  1 at Tbottom ¼2.2 1C). The resulting relationships are Zmax ¼ 3.956 (Basymp/N3)1/4 for the realizable k–ε turbulence model and Zmax ¼3.875 (Basymp/N3)1/4 for the RNG k–ε model.

rke

rngke

Fig. 10. Scaling of maximum plume rise height (Zmax) and neutrally buoyant plume height (Zneutral) with vent exit buoyancy flux (Bexit) and buoyancy frequency (N). (A) Determined from eight simulation runs (i.e. rke1-6, rke7_0 and rke7_1) listed in Table 2 using the realizable k–ε model (rke). (B) Determined from eight simulation runs (i.e. rngke1-6, rngke7_0 and rngke7_1) listed in Table 3 using the RNG k–ε model (rngke).

These relationships are similar to the commonly used scaling relationship (Eq. (1)). Thus, for applying correctly the scaling Eq. (1) in hydrothermal plume modeling, one needs to use the asymptotic buoyancy flux in Eq. (1). 3.4. Total plume vertically transported water mass flux scaling Using our numerical simulation data, we determine the vertical mass flux (M) inside the plume as a function of the vertical R1 position (z) according to MðzÞ  0 2πρwrdr (see the text following Eq. (3) for explanation of the variables). M increases (nonlinearly) from vent issued fluid mass flux (Mexit, Column 13 in

H. Jiang, J.A. Breier / Deep-Sea Research I 92 (2014) 41–55

51

Y=0.06765X+1.532 (R2=0.99873)

rke

Fig. 11. Scaling of total plume vertically transported water mass flux (Mtotal) with vent exit buoyancy flux (Bexit) and buoyancy frequency (N). Note that in the plot the vertical axis is Mtotal/ρbottom, where ρbottom is the density of the bottom seawater. The data are from eight simulation runs (i.e. rke1-6, rke7_0 and rke7_1) listed in Table 2 using the realizable k–ε model (rke).

Table 2) at the vent orifice level to its maximum at
0.71Zmax, and then decreases to zero at Zmax because of lateral spreading (figure not shown). The maximum M reached at
0.71Zmax vertical level is the total plume vertically transported water mass flux (Mtotal, Column 14 in Table 2) before the plume starts to spread laterally. We determine the scaling of Mtotal to be (Fig. 11): M total ¼ 0:068ρbottom ðBexit 3 N  5 Þ1=4 ;

ð8Þ

where ρbottom is the density of the bottom seawater. Entrainment and mixing caused by hydrothermal plumes is important for abyssal water mass transformation. The scaling relationship Eq. (8) may be used to quantify abyssal water mass flux conveyed upward by a hydrothermal plume. The scaling relationship confirms that the vertically transported abyssal water mass flux depends not only on hydrothermal flow injection at the vent orifice, but also more greatly on the background stratification (Columns 13 and 14 in Table 2). 3.5. Turbulent viscosity scaling Dimensional analysis of the simulation-derived turbulent viscosity data suggests a scaling for turbulent viscosity (νt) along the vent centerline (Fig. 12):
1=2 B νt ¼ 0:0170 exit f ðzn Þ ð9Þ N where zn is defined as the vertical distance above the vent orifice normalized by the maximum plume rise height (Zmax), and f(zn) is a function of zn that is shown by the solid line in Fig. 12A and tabulated by a text file (zstar_fzstar.txt) that is supplied as electronic supplementary material.

4. Discussion Our baseline model predictions appear to be consistent with the CTD observations of the neutrally buoyant plumes above the northern region of the ABE hydrothermal field and the rising plume in situ temperature measurements at 10 and 40 m above the ABE-A1 rising plume vent. The main discrepancy is the difference between the modeled and measured temperatures at the near vent 0.5 m elevation. As noted in Section 3.1.2, this discrepancy suggests that both/either (i) the temperature measurement at 0.5 m was closer to the periphery of the plume rather than the centerline, and/or (ii) some of our simplifying model assumptions and estimated parameters are not fully representative of the conditions at the time the measurements were made.

Y = 0.0170 X (R2 = 0.9714)

Fig. 12. Scaling of turbulent viscosity (νt) along the vent centerline with vent exit buoyancy flux (Bexit) and buoyancy frequency (N). (A) νt normalized by the maximum turbulent viscosity (νt,max) as a function of z* defined as (z  H)/Zmax. The triangle symbols are simulation grid point data from the eight simulation runs (i.e. rke1-6, rke7_0 and rke7_1) listed in Table 2 using the realizable k–ε model (rke), and the solid line is the averaged vertical profile of νt/νt,max. (B) Linear fitting to find the relationship between νt,max and (Bexit/N)1/2 using the data from the eight simulation runs shown in (A). The diamond symbols are simulation data points, and the solid line is the fit that is forced to pass through the (0, 0) point.

The former cannot be ruled out as plume motion and ROV position can be difficult to coordinate. However, the latter is also probable because our model greatly simplifies the actual environment, and this deserves discussion. The following are the explicit and implicit simplifying assumptions in our model and its comparison with the field observations. (1) The neutrally buoyant plumes identified in the CTD profiles are contiguous with the ABE-A1 plume. (2) Neighboring vent plumes are neglected. (3) Vent temperature, velocity, and salinity are constant in time. (4) Plume entrainment is constant in time. (5) Local bathymetry, other than the vent edifice, is neglected. (6) The vent is modeled as a single orifice. Regarding assumptions (1) and (2), there are at least 5 vent edifices within a 60 m radius of ABE-A1, a diffuse flow area on the elevated rise 100 m to the west, and many more vents to the south in the ABE field. Thus the CTD neutrally buoyant plume observations and the near vent ABE-A1 rising plume temperature measurements may not represent a contiguous plume in isolation. The mixing of multiple nearby plumes as they rise can result in an aggregate plume that may rise higher than the individual contributing plumes (Bemis et al., 2002; Rona et al., 2002). Regarding assumptions (3) and (4), vent temperature, velocity, and salinity and bottom currents, which can influence plume entrainment, are all variable (Rona et al., 2006). However, this variability is largely unconstrained in this field. Regarding assumption (5), ABE field topography is relatively complex. Specifically, ABE-A1 and nearby vent edifices sit upon finger like lava flows that provide up to 4 m of additional local bathymetric relief and a series of fault scarps create a
60 m E–W bathymetric gradient within the field (Ferrini et al., 2008). This bathymetric complexity is neglected in our model and may influence local bottom currents and potentially plume entrainment. Regarding assumption (6), modeling the vent as a single orifice rather than the actual multiple orifices may neglect some near vent mixing. Considering the limitations of these assumptions and the fact that both vent velocity and orifice area are estimates, a discrepancy between the model and observations is not unexpected. While the baseline model may be an imperfect replication of the ABE-A1 near vent observations, it does agree well with the other observations suggesting it is realistic for

H. Jiang, J.A. Breier / Deep-Sea Research I 92 (2014) 41–55

this area of the ABE field and thus is useful for understanding the physical structure of a hydrothermal plume. Our modeling shows how the steady-state hydrothermal plume flow pattern is characterized by a relatively constant entrainment coefficient (αe) until a plume elevation of
0.6Zmax (Fig. 7). Above that elevation, the most prominent flow feature is lateral spreading centered around
0.8Zmax. This is well described in previous research. Less well described are the two branches of inward flow above and below the lateral spreading flow. These inward flows partially compensate for the mass loss due to the strong lateral spreading flow that originates from the plume cap region. The lateral spreading flow consists of vent fluid and seawater entrained within the lower part of the rising plume, plus seawater delivered to the plume by the two branches of inward flow. Thus hydrothermal plumes are zones of mixing for water and material derived from the subseafloor, the benthic boundary layer, and the water column above and below the lateral spreading flow. The two branches of inward flow above and below the lateral spreading flow are predicted from numerical simulations. If they are real, they should reveal their existence in observational data and can be helpful in interpreting the observational data. Turbulent mixing is strengthened in the vicinity of hydrothermal plumes. First, turbulent mixing is strong in the region between
0.6Zmax and Zmax, as indicated by our model results for the distribution of the conservative scalar concentration (black line contours, Fig. 7A). This is consistent with the predictions that the largest values of turbulent viscosity (νt) and a strong recirculating flow occur in this same region (Fig. 6C and Fig. 7A). In addition, shear instabilities may also enhance mixing. These are likely to occur along the boundaries between the lateral spreading flow and the two branches of inward flow. Enhanced mixing both below and above the lateral spreading flow may greatly promote exchange of water mass and dispersal of plume chemical constituents between the lateral spreading plume and the background seawater. The dilution factor has been previously used to characterize the dilution of vent fluid by ambient seawater as the plume rises up before reaching the lateral spreading level (Rudnicki and Elderfield, 1992). For our baseline modeling case (rke7_0, Table 2), we calculate the dilution factor as the ratio between the vent orifice issued fluid mass flux (Mexit, Column 13 in Table 2) and the total plume vertically transported water mass flux (Mtotal, Column 14 in Table 2). The resulting dilution factor is
1.43  10  4, indicating that the vent fluid is diluted
7000 times with ambient seawater at
0.71Zmax vertical level beyond which the lateral spreading occurs. However, the concentration of a vent orifice issued conservative scalar is a more accurate indicator of the dilution (black line contours, Fig. 7A). At
0.71Zmax and within a much narrower region than the plume velocity width at that height, the concentration of the conservative scalar remains 4 7  10  3, which is much larger than the predicted dilution factor. In other words, the scalar width is much thinner than the plume velocity width. Thus, the dilution factor is probably too “averaged out” to represent the actual spatial distribution and heterogeneity of the dilution of vent fluid by ambient seawater. The k–ε turbulence closure model calculates turbulent scalar diffusivity (KC) according to KC ¼ νt/Sct, where Sct is the turbulent Schmidt number ( ¼0.7). KC is calculated in the range of
7.2  10  3 to 9.8  10  2 m2 s  1 in the simulated plume cap region. This is similar to previous results. Walter et al. (2010) estimated an averaged turbulent diffusivity of 4  10  2 m2 s  1 close to the vent site of the hydrothermal plume observed at the Nibelungen field (81180 S 131300 W) by analyzing temperature finestructure data. Keir et al. (2008) calculated a vertical turbulent diffusivity in the range of 9  10  3–8  10  2 m2 s  1 in order to

reproduce the observed distribution of gases emitted from the Drachenschlund hydrothermal vent (81180 S 131300 W) using a Gaussian plume model. Our numerical simulations suggest a scaling relationship for the maximum turbulent scalar diffusivity along the vent centerline with respect to Bexit and N:
K C; max ¼ 0:024

Bexit N

1=2 ;

ð10Þ

which is based on Eq. (9) with an assumption of Sct ¼0.7. KC,max scales with (Bexit/N)1/2 while Zmax scales with (Bexit/N3)1/4. These two differential scaling relationships control the steady state concentration of an inert species issued from the vent orifice (e.g. helium-3) within the plume. Specifically, the steady state concentration falls more rapidly along the vent centerline (normalized by Zmax) with decreasing either Bexit or N (Fig. 13). Thus, the vent condition (Bexit) and the background stratification (N) jointly put a predicable physical control on the vertical distribution of vent derived chemical species in the plume. Our simulation-derived steady-state turbulence dissipation rate (ε) field attains its maximum value close to the vent orifice (Fig. 6B). By contrast, our simulation-derived steady-state turbulent viscosity (νt) field reaches its maximum value close to the maximum plume rise height (Fig. 6C). Based on a simple balance between local shear production, dissipation, and loss to buoyancy derived from the turbulent kinetic energy equation, Osborn (1980) proposed to use Kρ r0.2  ε/N2 to estimate turbulent diffusivity from turbulence dissipation measurements. If the Osborn model were used to calculate turbulent diffusivity from our simulationderived turbulence dissipation rate, it would generate a turbulent diffusivity field that is much different from that calculated from using KC ¼ νt/Sct. We therefore suggest that the Osborn model may not be appropriate for hydrothermal plumes where buoyancy contributes significantly to the production of turbulent kinetic energy. Significant buoyancy contribution to turbulent kinetic

rrke7_1 ke7 _3

52

1.2

0 7_ ke6 e3 e1 r rk rk

rke

1.0

1.8

0.6

0.4

0.2

0.0 0.001

0.010

0.100

1.000

Fig. 13. Steady state scalar concentration, C, profiles along the vent centerline, plotted against the normalized vertical distance above the vent exit, (z  H)/Zmax. C¼1 at the vent orifice and C¼0 at the top boundary. The C profiles are shown to vary with background buoyancy frequency (N; rke1, rke3, rke6 and rke7_0 of Table 2), vent exit buoyancy flux (Bexit; rke7_0 and rke7_1 of Table 2), and vent chimney height above seafloor (H; rke7_0 and rke7_3 of Table 2).

H. Jiang, J.A. Breier / Deep-Sea Research I 92 (2014) 41–55

energy has implications for the spatial and temporal characteristics of turbulent mixing and diffusivity. Turbulent viscosity (νt) can be written as the product of a turbulent velocity scale (w0 ) and a turbulent length scale (l0 ), i.e. νt ¼ w0  l0 . Because the steady-state turbulent kinetic energy (k) decreases with increasing height above the vent orifice (Fig. 6A) and because w0 p k1/2, w0 must also decrease with height. On the other hand, νt increases with height (Fig. 6C). Thus, the turbulent length scale (l0 ) must increase with height in a greater pace than that in which w0 decreases with height. This is phenomenologically consistent with hydrothermal plume observation that the size of smoke billows increases with height. Turbulent entrainment of ambient seawater into the rising plume plays a key role in the hydrodynamics of the hydrothermal plume. It is thus vital that the turbulence closure model used in our simulations generates the right amount of entrainment, quantified by the entrainment coefficient (αe). Experimental studies have suggested αe to be in the range of 0.10–0.16 for buoyancy-driven plumes (Fischer et al., 1979; Chen and Rodi, 1980; Wang and Law, 2002). Specifically in the case of hydrothermal plumes, in situ acoustic imaging of the lower 25 m section of a buoyant hydrothermal plume, has been used to estimate that αe was in the range of 0.07–0.18 (Rona et al., 2006). The highest values of αe in the Rona et al. (2006) study was attributed to the case where the plume was “bent” by background flow, which was not a case considered by the present study. In our simulations, the realizable k–ε model (rke) generates αe in the range of 0.10–0.16 (maximum standard deviation 0.020) below the vertical level of
0.6Zmax (Fig. 8A). By contrast, the RNG k–ε model (rngke) generates larger αe values in the range of 0.10–0.19 (maximum standard deviation 0.025, Fig. 8B). Thus based on the available data both of these models produce reasonable values for αe, but the realizable k–ε model (rke) produces a more specific match to observations. On this basis we have designated simulation rke7_0 (Table 2, using the realizable k–ε model) as our baseline modeling case. The realizable k–ε model may actually perform better than the RNG k–ε model for modeling the hydrothermal plume flow. It is known that the RNG k–ε model overpredicts the spreading rate of the round jet. By contrast, the realizable k–ε model has been shown to predict round jet spreading correctly (Shih et al., 1995). Our simulations indicate that vent fluid exit velocity and vent chimney height both affect the local entrainment in the near vent region (Fig. 9A and B). The affected region is
10 times the orifice diameter surrounding the vent. Many of the initial chemical reactions that occur when high temperature vent fluid, supersaturated in inorganic constituents, is quenched by seawater occur in this zone. As these reactions are by their nature sensitive to both temperature and mixing, the rate and extent to which they proceed will also be affected by near vent physical conditions. As these initial reactions establish a starting point for subsequent hydrothermal plume biotic and abiotic processes, it follows that near vent physical conditions can influence the biogeochemistry of the rising plume as well. These simulation results are being used to drive models of coupled microbial-geochemical dynamics in deep-sea hydrothermal plumes (Breier et al., 2011).

5. Conclusions The present study shows that it is feasible to use a CFD approach to generate a numerical simulation that matches measured vent conditions (such as vent geometry, exit fluid temperature and velocity, and background stratification) and reproduces a maximum plume rise height estimated from observation. Such a simulation provides data for properties that cannot be measured directly, i.e. the plume flow velocity field, temperature field, and

53

turbulent viscosity field. These numerical data are particularly detailed for the near vent region where important plume chemical and biogeochemical processes are likely to occur. Thus, this CFD model can be used to provide the physical driving forces for the development of numerical modeling of chemical and biological processes in a hydrothermal plume. The present high-resolution CFD simulation shows that the vertical distribution of mixing and entrainment in the plume is more complex than predicted by analytical theories. The classical entrainment assumptions are shown to only apply up to the vertical level of
0.6 times the maximum plume rise height (Zmax). Below that level, the entrainment coefficient remains relatively constant (
0.15). Above that level, the plume flow consists of a pronounced lateral spreading flow, two branches of inward flow immediately above and below the lateral spreading, and recirculation flanking the plume cap region. Both turbulent kinetic energy and turbulence dissipation rate reach their maximum near the vent; however, turbulent viscosity attains its maximum near the plume top, indicating strong turbulent mixing in that region. Near vent physical conditions, including chimney height and fluid exit velocity, are also shown to influence plume mixing from the vent orifice to a distance of
10 times the vent orifice diameter. Thus, physical parameters place a strong kinetic constraint on the chemical reactions occurring in the initial particle-forming zone of a hydrothermal plume. Although CFD simulations enable mixing and transport predictions with fewer simplifying assumptions, reliable predictions may be hampered by the present difficulty in obtaining measurements of near vent parameters such as orifice area and fluid velocity. The CFD simulation approach allows for deriving several important scalings based on background stratification (N) and the actual plume exit buoyancy flux (Bexit) that is calculated directly from the thermal and kinematic properties of the exit vent fluid. Maximum plume rise height (Zmax) scales as Zmax ¼2.0 (Bexit/N3)1/4. Neutrally buoyant plume height (Zneutral) scales as Zneutral ¼1.5 (Bexit/N3)1/4. Maximum turbulent diffusivity (KC,max) scales as KC,max ¼ 0.024(Bexit/N)1/2. Total plume vertically transported water mass flux (Mtotal, i.e. before the plume starts to spread laterally) scales as Mtotal ¼0.068  ρbottom(B3exit/N5)1/4, where ρbottom is the density of the bottom seawater. These simple scaling relationships relate plume entrainment, mixing and transport parameters to their actual physical controls. These scalings can be useful for constructing 1-dimensional models for geochemistry and microbial activity in a hydrothermal plume. The scaling relationships of maximum plume rise height and total plume vertically transported water mass flux can be useful for understanding the role of deep-sea hydrothermal plumes in transporting and mixing abyssal water mass. The results of the present numerical simulation highlight some important features of hydrothermal plume hydrodynamics that are especially relevant to plume impact on the ocean.

Acknowledgements Thanks to C. Fisher, P. Girguis, G. Luther III, A. L. Reysenbach, M. Tivey, and the Eastern Lau Spreading Center 2009 scientific parties (NSF OCE-0424953, OCE02-040985, OCE-0728391, OCE-0752469, OCE-0751839) for allowing JAB to participate in their cruises. Thanks also to G. Luther III, A. Madison, J. Sylvan, K. Edwards, and S. White for collaborating with JAB during this field work and providing measurements used in this study. This work was supported by a National Science Foundation grant NSF OCE1038055 through the RIDGE2000 program. This work was also funded by the Gordon and Betty Moore Foundation through Grant GBMF2609 to Dr. G. J. Dick at University of Michigan (subcontract

54

H. Jiang, J.A. Breier / Deep-Sea Research I 92 (2014) 41–55

to JAB and HJ). We thank two anonymous reviewers for providing very helpful and constructive comments that improved the manuscript. We thank G. J. Dick, D. C. Reed and B. M. Toner for very helpful discussions. HJ is grateful to J. Lin for very stimulating discussions.

Appendix A. Supporting information Supplementary data associated with this article can be found in the online version at http://dx.doi.org/10.1016/j.dsr.2014.06.006. References Baker, E.T., 2007. Hydrothermal cooling of midocean ridge axes: do measured and modeled heat fluxes agree? Earth Planet. Sci. Lett 263, 140–150. Baker, E.T., German, C.R., Elderfield, H., 1995. Hydrothermal plumes over spreadingcenter axes: global distributions and geological inferences. In: Humphris, S.E., Zierenberg, R.A., Mullineaux, L.S., Thomson, R.E. (Eds.), Seafloor Hydrothermal Systems: Physical, Chemical, Biological, and Geological Interactions. Geophysical Monograph, 91. American Geophysical Union, pp. 47–71. Baker, E.T., Lavelle, J.W., Feely, R.A., Massoth, G.J., Walker, S.L., Lupton, J.E., 1989. Episodic venting of hydrothermal fluids from the Juan de Fuca Ridge. J. Geophys. Res. 94 (B7), 9237–9250. Bemis, K., Rona, P., Jackson, D., Jones, C., Silver, D., Mitsuzawa, K., 2002. A comparison of black smoker hydrothermal plume behavior at Monolith Vent and at Clam Acres Vent field: dependence on source configuration. Mar. Geophys. Res. 23, 81–96. Bennett, S.A., Achterberg, E.P., Connelly, D.P., Statham, P.J., Fones, G.R., German, C.R., 2008. The distribution and stabilisation of dissolved Fe in deep-sea hydrothermal plumes. Earth Planet. Sci. Lett. 270, 157–167. Breier, J.A., Toner, B.M., Fakra, S.C., Marcus, M.A., White, S.N., Thurnherr, A.M., German, C.R., 2012. Sulfur, sulfides, oxides and organic matter aggregated in submarine hydrothermal plumes at 91500 N East Pacific Rise. Geochim. Cosmochim. Acta 88, 216–236. Breier, J.A., Osicki, O., Jiang, H., Anantharaman, K., Dick, G., Wendt, K., Sorenson, J.V., Toner, B., 2011. Mineral formation and trace element uptake in rising hydrothermal plumes of the Lau basin, In: Abstract OS23B-01 Presented at 2011 Fall Meeting, AGU, San Francisco, CA., 5–9, Dec. Briggs, G.A., 1969. Optimum formulas for buoyant plume rise. Philos. Trans. R. Soc. London, Ser. A 265, 197–203. Cann, J.R., Strens, M.R., 1989. Modeling periodic megaplume emission by black smoker systems. J. Geophys. Res. 94 (B9), 12227–12237. Chen, J.C., Rodi, W., 1980. Vertical Turbulent Buoyant Jets—A Review of Experimental Data. Pergamon Press, New York, NY. Converse, D.R., Holland, H.D., Edmond, J.M., 1984. Flow rates in the axial hot springs of the East Pacific Rise (211N): implications for the heat budget and the formation of massive sulfide deposits. Earth Planet. Sci. Lett. 69, 159–175. Cowen, J.P., Massoth, G.J., Baker, E.T., 1986. Bacterial scavenging of Mn and Fe in a mid-field to far-field hydrothermal particle plume. Nature 322, 169–171. de Angelis, M.A., Lilley, M.D., Baross, J.A., 1993. Methane oxidation in deep-sea hydrothermal plumes of the Endeavour segment of the Juan de Fuca Ridge. Deep Sea Res. Part I, 1169–1186 (Papers 40). Dymond, J., Roth, S., 1988. Plume dispersed hydrothermal particles: a time-series record of settling flux from the Endeavour Ridge using moored sensors. Geochim. Cosmochim. Acta 52, 2525–2536. Elderfield, H., Schultz, A., 1996. Mid-ocean ridge hydrothermal fluxes and the chemical composition of the ocean. Annu. Rev. Earth Planet. Sci. 24, 191–224. Ferrini, V.L., Tivey, M.K., Carbotte, S.M., Martinez, F., Roman, C., 2008. Variable morphologic expression of volcanic, tectonic, and hydrothermal processes at six hydrothermal vent fields in the Lau back-arc basin. Geochem. Geophys. Geosyst. 9, Q07022, http://dx.doi.org/10.1029/2008GC002047. Fischer, H.B., List, E.J., Imberger, J., Brooks, N.H., 1979. Mixing in Inland and Coastal Waters. Academic, San Diego. Fofonoff, N.P., Millard Jr., R.C., 1983. Algorithms for computation of fundamental properties of seawater. UNESCO Tech. Pap. Mar. Sci., 44. Kadko, D., 1993. An assessment of the effect of chemical scavenging within submarine hydrothermal plumes upon ocean geochemistry. Earth Planet. Sci. Lett. 120, 361–374. Kim, S.L., Mullineaux, L.S., Helfrich, K.R., 1994. Larval dispersal via entrainment into hydrothermal vent plumes. J. Geophys. Res. 99 (C6), 12655–12665. Keir, R.S., Schmale, O., Walter, M., Sültenfuß, J., Seifert, R., Rhein, M., 2008. Flux and dispersion of gases from the “Drachenschlund” hydrothermal vent at 81180 S, 131300 W on the Mid-Atlantic Ridge. Earth Planet. Sci. Lett. 270, 338–348. Lavelle, J.W., 1997. Buoyancy-driven plumes in rotating, stratified cross flows: plume dependence on rotation, turbulent mixing, and cross-flow strength. J. Geophys. Res. 102 (C2), 3405–3420. Lavelle, J.W., 1995. The initial rise of a hydrothermal plume from a line segment source—results from a three-dimensional numerical model. Geophys. Res. Lett. 22, 159–162. Lavelle, J.W., 1994. A Convection Model for Hydrothermal Plumes in a Cross Flow. NOAA Technical Memorandum ERL PMEL-102, NTIS: PB94-218815, 18.

Lavelle, J.W., Baker, E.T., 1994. A numerical study of local convection in the benthic ocean induced by episodic hydrothermal discharges. J. Geophys. Res. 99 (C8), 16065–16080. Lavelle, J.W., Di Iorio, D., Rona, P., 2013. A turbulent convection model with an observational context for a deep-sea hydrothermal plume in a time-variable cross flow. J. Geophys. Res. Oceans 118, 1–16. Lesniewski, R.A., Jain, S., Schloss, P.D., Anantharaman, K., Dick, G.J., 2012. The metatranscriptome of a deep-sea hydrothermal plume is dominated by water column methanotrophs and chemolithotrophs. ISME J. 6, 2257–2268. Li, M., Toner, B.M., Baker, B.J., Breier, J.A., Sheik, C.S., Dick, G.J., 2014. Microbial iron uptake as a mechanism for dispersing iron from deep-sea hydrothermal vents. Nat. Commun., 5, http://dx.doi.org/10.1038/ncomms4192. Little, S.A., Stolzenbach, K.D., Von Herzen, R.P., 1987. Measurements of plume flow from a hydrothermal vent field. J. Geophys. Res. 92 (B3), 2587–2596. Lupton, J.E., Delaney, J.R., Johnson, H.P., Tivey, M.K., 1985. Entrainment and vertical transport of deep-ocean water by buoyant hydrothermal plumes. Nature 316, 621–623. Luther III, G.W., Glazer, B.T., Ma, S., Trouwborst, R.E., Moore, T.S., Metzger, E., Kraiya, C., Waite, T.J., Druschel, G., Sundby, B.R., Taillefert, M., Nuzzio, D.B., Shank, T.M., Lewis, B.L., Brendel, P.J., 2008. Use of voltammetric solid-state (micro)electrodes for studying biogeochemical processes: laboratory measurements to real time measurements with an in situ electrochemical analyzer (ISEA). Mar. Chem. 108, 221–235. McDougall, T.J., 1990. Bulk properties of “hot smoker” plumes. Earth Planet. Sci. Lett. 99, 185–194. McDuff, R.E., 1995. Physical dynamics of deep-sea hydrothermal plumes. In: Humphris, S.E., Zierenberg, R.A., Mullineaux, L.S., Thomson, R.E. (Eds.), Seafloor Hydrothermal Systems: Physical, Chemical, Biological, and Geological Interactions. Geophysical Monograph, 91. American Geophysical Union, Washington, D.C., pp. 357–368. Middleton, J.M., Thomson, R.E., 1986. Modelling the rise of hydrothermal plumes. Can. Tech. Rep. Hydrogr. Ocean Sci. 69, 18. Morton, B.R., Taylor, G.I., Turner, J.S., 1956. Turbulent gravitational convection from maintained and instantaneous sources. Proc. R. Soc. London, Ser. A 234, 1–23. Mottl, M.J., Seewald, J.S., Wheat, C.G., Tivey, M.K., Michael, P.J., Proskurowski, G., McCollom, T.M., Reeves, E., Sharkey, J., You, C.-F., Chan, L.-H., Pichler, T., 2011. Chemistry of hot springs along the Eastern Lau Spreading Center. Geochim. Cosmochim. Acta 75, 1013–1038. Osborn, T., 1980. Estimates of the local rate of vertical diffusion from dissipation measurements. J. Phys. Oceanogr. 10, 83–89. Pham, M.V., Plourde, F., Kim, S.D., 2006. Large-eddy simulation of a pure thermal plume under rotating conditions. Phys. Fluids 18, 015101, http://dx.doi.org/ 10.1063/1.2162186. Plourde, F., Pham, M.V., Kim, S.D., Balachandar, S., 2008. Direct numerical simulations of a rapidly expanding thermal plume: structure and entrainment interaction. J. Fluid Mech. 604, 99–123. Ramondenc, P., Germanovich, L.N., Von Damm, K.L., Lowell, R.P., 2006. The first measurements of hydrothermal heat output at 91500 N, East Pacific Rise. Earth Planet. Sci. Lett. 245, 487–497. Rona, P.A., Bemis, K.G., Jones, C.D., Jackson, D.R., Mitsuzawa, K., Silver, D., 2006. Entrainment and bending in a major hydrothermal plume, Main Endeavour Field, Juan de Fuca Ridge. Geophys. Res. Lett. 33, L19313, http://dx.doi.org/ 10.1029/2006GL027211. Rona, P.A., Bemis, K.G., Silver, D., Jones, C.D., 2002. Acoustic imaging, visualization, and quantification of buoyant hydrothermal plumes in the ocean. Mar. Geophys. Res. 23, 147–168. Rudnicki, M.D., Elderfield, H., 1992. Theory applied to the Mid-Atlantic Ridge hydrothermal plumes: the finite-difference approach. J. Volcanol. Geotherm. Res. 50, 161–172. Sander, S.G., Koschinsky, A., Massoth, G., Stott, M., Hunter, K.A., 2007. Organic complexation of copper in deep-sea hydrothermal vent systems. Environ. Chem. 4, 81–89. Shih, T.-H., Liou, W.W., Shabbir, A., Yang, Z., Zhu, J., 1995. A new eddy-viscosity model for high Reynolds number turbulent flows—model development and validation. Comput. Fluids 24, 227–238. Speer, K.G., Marshall, J., 1995. The growth of convective plumes at seafloor hot springs. J. Mar. Res. 53, 1025–1057. Speer, K.G., Rona, P.A., 1989. A model of an Atlantic and Pacific hydrothermal plume. J. Geophys. Res. 94 (C5), 6213–6220. Stranne, C., Sohn, R.A., Liljebladh, B., Nakamura, K., 2010. Analysis and modeling of hydrothermal plume data acquired from the 851E segment of the Gakkel Ridge. J. Geophys. Res. 115, C06028, http://dx.doi.org/10.1029/2009JC005776. Sun, H., Feistel, R., Koch, M., Markoe, A., 2008. New equations for density, entropy, heat capacity, and potential temperature of a saline thermal fluid. Deep Sea Res. Part I 55, 1304–1310. Suzuki, Y.J., Koyaguchi, T., 2010. Numerical determination of the efficiency of entrainment in volcanic eruption columns. Geophys. Res. Lett. 37, L05302, http://dx.doi.org/10.1029/2009GL042159. Tagliabue, A., Bopp, L., Dutay, J.-C., Bowie, A.R., Chever, F., Jean-Baptiste, P., Bucciarelli, E., Lannuzel, D., Remenyi, T., Sarthou, G., Aumont, O., Gehlen, M., Jeandel, C., 2010. Hydrothermal contribution to the oceanic dissolved iron inventory. Nat. Geosci. 3, 252–256. Tao, Y., Rosswog, S., Brüggen, M., 2013. A simulation modeling approach to hydrothermal plumes and its comparison to analytical models. Ocean Modell. 61, 68–80. Thorpe, S.A., 2005. The Turbulent Ocean. Cambridge University Press, Cambridge.

H. Jiang, J.A. Breier / Deep-Sea Research I 92 (2014) 41–55

Thurnherr, A.M., St. Laurent, L.C., 2012. Turbulence observations in a buoyant hydrothermal plume on the East Pacific Rise. Oceanography 25, 180–181. Toner, B.M., Fakra, S.C., Manganini, S.J., Santelli, C.M., Marcus, M.A., Moffett, J.W., Rouxel, O., German, C.R., Edwards, K.J., 2009. Preservation of iron(II) by carbonrich matrices in a hydrothermal plume. Nat. Geosci. 2, 197–201. Turner, J.S., 1973. Buoyancy Effects in Fluids. Cambridge University Press, Cambridge. Turner, J.S., Campbell, I.H., 1987. Temperature, density and buoyant fluxes in “black smoker” plumes, and the criterion for buoyancy reversal. Earth Planet. Sci. Lett. 86, 85–92.

55

Walter, M., Mertens, C., Stöber, U., German, C.R., Yoerger, D.R., Sültenfuß, J., Rhein, M., Melchert, B., Baker, E.T., 2010. Rapid dispersal of a hydrothermal plume by turbulent mixing. Deep Sea Res. Part I 57, 931–945. Wang, H., Law, A.W.-K., 2002. Second-order integral model for a round turbulent buoyant jet. J. Fluid Mech. 459, 397–428. Woods, A.W., 2010. Turbulent plumes in nature. Annu. Rev. Fluid Mech. 42, 391–412.