The model describes a thermal plume rising into a temperature- and salinity-stratified background environment marked by steady, vertically sheared cross flow. The fluid is incompressible and the response to the localized buoyancy source nonhydrostatic. Together with the equation of state for seawater [Fofonoff and Millard, 1983], the underlying model equations in the Boussinesq approximation are
(1) |
(2) |
(3) |
(4) |
(5) |
(6) |
(7) |
(8) |
(9) |
(10) |
where t is time; consists of u, v, and vertical w (positive upward) components; p is pressure; is density; (1028.11 kg m) is a reference density; g (9.81 m s) is acceleration of gravity; is the rotation vector; and is the unit vector in the vertical direction. Potential temperature, , and salinity, S, depend on rates of discharge of heat QH, via Q = QH /(CP), and salt QS, where CP (4200 J kg °C) is the specific heat of seawater. Equations (1)-(3) are statements of conservation for momentum, mass (given fluid incompressibility), and heat (or salt). Chemical tracers have conservation equations analogous to (3). Retention of all terms in the vertical force balance in (1) makes the model nonhydrostatic. The viscosity tensor A and diffusivity tensor K are time- and space-dependent; diffusive transport is down gradient.
Equations (4)-(10) constitute the turbulent mixing submodel. Horizontal and vertical components of mixing tensors AH or V (and KH or V via (10)) are each composed of two parts: an isotropic (AI) mixing coefficient and either horizontal (AHMIN) or vertical (AVMIN) background terms (equations (4)-(5)). The isotropic mixing coefficient AI (equation (6)) is made to depend on fluid shear (Sij, here given in Cartesian form, (9)), on shear Richardson number (RI (8)), on Prandtl number (PR (10)), and on the turbulence length scale (lS) and Smagorinsky (CS) constants (equation (6)). If RI were zero, (6) would reduce to a form originating with Smagorinsky [1963]; it stems from an assumed local balance between shear production and turbulence dissipation in the turbulence kinetic energy (TKE) equation. When buoyancy production or dissipation is added to the local TKE balance, the RI dependent factor becomes part of the formulation [Lilly, 1962]. Advection and diffusion of TKE are generally found to represent small contributions to the TKE equation balance when explicitly evaluated, and they are here implicitly set to zero. Mixing terms have no explicit dependence on rotation. The reason is that in forming TKE, terms involving sum to zero.
Smagorinsky-Lilly mixing, with or without the RI dependent factor, is in widespread use in atmospheric calculations [e.g., Clark and Farley, 1984; Mason, 1989]. Subgrid-scale mixing of this type is frequently used in large eddy simulations (LES) of both engineering and geophysical natures [e.g., Reynolds, 1990; Galperin and Orszag; 1993; Smagorinsky, 1993; Mason, 1994]. Following the nomenclature of Reynolds [1990], Wyngaard [1990] notes that models using this kind of turbulence closure should actually be typed very large eddy simulation (VLES). VLES implies that, although resolved scales include much of the eddy energy, model resolution is not high enough that smallest resolved scales fall within the subinertial range of turbulence where details of dissipation ought not to matter.
Background mixing terms (equations (4)-(5)), presumed to be small, have here been added to AI to ensure a modest amount of mixing in regions without shear, e.g., away from the plume region, where AI will be zero. Other closure models are in current use as well. Sykes et al. [1986] used a full TKE equation to determine q, the square root of twice TKE; turbulent mixing was then made proportional to the product of q and lS. In both the present, Sykes et al.'s, and many other approaches, lS has a size that is imposed rather than calculated. The model summarized by (1)-(10), but without cross flow, has been studied in the context of laboratory convection experiments as a way of investigating the value of the subgrid turbulence closure subcomponent [Lavelle and Smith, 1996].
For hydrothermal venting, convective flow is forced by a source of heat at interior points along the lower boundary of the region of interest. Steady background cross flow is forced by a constant horizontal pressure gradient, P. Without point source heating, the steady background velocity profile uBKG = (uBKG, vBKG) results from an Ekman horizontal force balance (equation (1)):
- P/ - 2 × uBKG + d/dz [Az d(uBKG)/dz] = 0 | (11) |
and boundary conditions uBKG = 0 (z = 0, bottom) and duBKG/dz = 0 (z = h, top). With uniform P in the y direction, nonzero , and Az profiles that increase toward the seafloor as described later, flow is oriented in the x direction except near the bed. Profiles uBKG and vBKG (e.g., Figure 1a) were used to initialize velocity fields and maintain upstream velocity boundary conditions thereafter. Without source heat (and/or salt), uBKG and vBKG are maintained numerically over time as solutions to (1)-(10) as they must be. In reference experiment 21 (Table 1), outside the boundary layer, uBKG was 1.5 cm s, while vBKG was nonzero only within that layer (Figure 1a). Mean currents of that magnitude are typical at ridge crest depths on the JDFR [e.g., Cannon et al., 1991].
Figure 1. Idealized environmental profiles for u, v, AZ, , S, and . Profile shapes and magnitudes nominally represent conditions in the lowest 300 m of water column above the seafloor spreading center of the Juan de Fuca Ridge, northeast Pacific. The AZ profile controlling the boundary layer thickness is hypothetical, but the shape and magnitude are generally consistent with AZ values deduced from benthic ocean measurements at other sites.
Table 1. Parameter Values for the Convective Plume Experiments
Exp. No. | 2Z, s |
CS | AH
MIN, cm s |
U, m s |
h,
m |
21 | 1.03 × 10 | 0.2 | 100 | 0.015 | 300 |
22 | 0 | 0.2 | 100 | 0.015 | 300 |
23 | 1.03 × 10 | 0 | 100 | 0.015 | 300 |
26 | 1.03 × 10 | 0.2 | 10 | 0.015 | 300 |
24 | 1.03 × 10 | 0.2 | 100 | 0.03 | 200 |
25 | 1.03 × 10 | 0.2 | 100 | 0.06 | 150 |
The vertical force balance without heating is hydrostatic: dP/dz = gBKG, where BKG is background density; wBKG = 0. The variable p in (1) thus represents the sum of P(z), P(x,y), and a perturbation pressure p(x,y,z) caused by heating. In actual practice, hydrostatic terms are first subtracted from (1), with the consequence that in resulting equations, p represents the sum of P and p and becomes = ( - BKG ).
Background profiles (Figures 1a-1b) hinge on the form given Az. To simplify the initialization process and make background profiles independent of boundary layer shear, of (4) was set to zero. Vertical turbulent mixing coefficients are thus fixed from the start of calculations. Because vertical mixing rates in the ocean are small in the interior (AZMIN) and increase (to AZMAX) as z 0 [e.g., Garrett, 1979], the profile for A (Figure 1a) was given the form
(12) |
where hB determines the boundary layer thickness. The value of Az declines to within 1% of AZMIN by a distance of 5 hB of the seafloor (z = 0). Values of AZMIN = 1 cm s, AZMAX = 100 cm s, and a viscosity height scale length hB = 10 m were taken as nominally representative of benthic ocean conditions.
Without heating, background profiles BKG and SBKG must be time-invariant and laterally uniform. Since wBKG = 0, (3) in that case leads to
d/dz [K (dBKG/dz)] = 0 | (13) |
or once integrated over depth interval [z,h],
K [dBKG/dz] = F0 | (14) |
where F0 is a constant the value of which is set by conditions at z = h. Over depth interval [0,h], profiles of Kz and BKG are thus linked; where K grows large near the seafloor, dBKG/dz must adjust appropriately. Equations (13)-(14) imply that to maintain a steady BKG profile, the gradient flux of BKG must be balanced by a vertical flux of opposite sign. Otherwise, background profiles would continually change, making more difficult the resolution of small temperature anomalies caused only by hydrothermal heating. Analogous arguments apply to SBKG.
Above a boundary layer, background profiles of and S were taken to be linear: = (1 + ) and S = S0 (1 + ) , where , S, , and were given the values 1.6339°C, 34.616, 1.8925 × 10 m, and -1.5751 × 10 m and where is height above the sea floor. These values are taken from linear least squares fits of hydrographic profile data from the JDFR over the 2100-2400 m depth range, profiles judged least affected by hydrothermal influences. Given those profiles and (14), F0 = 0K (z = h) and F0S = 2SK (z = h). Since A(z = h) = 10 m s (equation (12)) and PR = 1.0 (equation (10)), F0 = 3.1 × 10 °C m s and F0S = 5.5 × 10 m s.
Integrating (4) with the A profile of (12) for PR = 1.0 results in the background profile BKG for 0 < z < h:
(15) |
having a well-mixed boundary layer (Figure 1b), the height (~50 m) of which is determined by hB. In (15), K represents (KZMAX - KZMIN). A like equation can be written for SBKG. For z >> hB, the second term of (15) is negligible. As with velocity, BKG and SBKG profiles (Figure 1b) are used to initialize fields and provide upstream boundary profiles. Note in Figure 1 that where gradients of and S are large, velocity shear is small, and vice versa. This situation suppresses shear instabilities that might otherwise develop. It follows that = - 1000 has a similarly well-mixed boundary layer (Figure 1b). Buoyancy frequency, N, over depths 2100-2350 m is 7.9 × 10 s.
A resolved Ekman boundary layer serves these purposes. It allows for directional shear of flow that can influence plume dispersion. It reduces flow speed in the boundary layer, flow that can affect entrainment of fluid into the plume stem, particularly on the downstream side. Finally, a resolved boundary layer is particularly important for plumes the B of which is not great enough to cause the plumes to rise distinctly above the seafloor. Diffuse low-temperature sources of hydrothermal heat are instances of that situation [Trivett and Williams, 1994].
Equations (1)-(10) were integrated over a rectangular domain with boundaries open in x, cyclic in y, and fixed in z. Specifically, conditions were
Lateral inflow boundary, x = X1
u = uBKG v = vBKG w = 0 = BKG S = SBKG. |
(16) |
Lateral outflow boundary, x = X
uXXX = vXX = wXX = X = SX = 0. | (17) |
Lower boundary, z = 0
u = v = w = 0 Kz [/z] = F0 Kz [S/z] = F0S. |
(18) |
Upper boundary, z = h
u/z
= v/z
= w = 0 Kz [/z] = F0 Kz [S/z] = F0S. |
(19) |
Front and back boundaries, y = Y and y = Y
(20) |
Boundary conditions on p in the x and z directions were taken in the manner of Harlow and Welch [1965]. Cyclic conditions were used in the y direction for all variables (u, v, w, , S, and p) for two reasons: so the Ekman boundary layer of the background flow would be uniform across the region of calculation and so plumes could pass though y boundaries as if transparent. Cross-flow strength and width of the solution domain were always large enough that the downstream plume never filled the domain side to side, though in some cases the plume passed through y boundaries. Boundary openness in the x direction was essential to allow passage of heat and tracers out of the domain, which if prevented would disallow any thermal equilibrium state from being reached. X direction boundary conditions are those suggested by Johansson [1993]. In these experiments, no absorbing layers [e.g., Clark and Farley, 1984] were used to damp waves.
As with background profiles, length scales, B, and background were chosen with regard to conditions for chronic hydrothermal plumes at the JDFR. Equations were solved in a Cartesian domain 640 × 320 × 300 m when U was 1.5 cm s, but domain height was reduced to 200 and then 150 m for increasingly larger U0 (Table 1). Domain depth h was chosen for each U so that plumes penetrated to ~ 0.5 h. x, y, z, t were 5, 5, and 9.4 (6.2 or 4.7) m, and 15 s, respectively. N, Ny, Nz, and N were 128, 64, 32, and 2000 (or 5760).
The source region was 10 × 10 m and was centered at the horizontal location (x,y) = (-105, 0) m; source lateral dimension will hereafter be represented by D. Individual hydrothermal vents have smaller surface cross sections, but separate vents can be clustered in fields [Ginster et al., 1994], or the hydrothermal source may be diffuse, like the sulfide mounds described by Schultz et al. [1992]. The intention here was to model plumes not from a single vent, but from a small composite venting source or vent field. Resolving a single small vent and its regional plume is beyond present computational resources. Model source area was thus chosen to be representative of typical vent field and sulfide mound surface areas. Despite having an extended area, the source is nonetheless point-like in that rise height, hRISE, is very much greater than D.
Source heat flux, Q, was set at 1.3 × 10 J s (13 MW), a value chosen to cause maximum plume rise of ~200 m under reference experiment conditions. Since buoyancy flux B = gQ//CP, where is the thermal expansion coefficient (7.3 × 10°C), B was fixed at 2.1 × 10 m s. Buoyancy flux density was consequently 2.1 × 10 m s.
The value assigned QH is consistent with the wide range of heat fluxes measured for single vents when extrapolated to a small vent field. For 18 vents on the southern Juan de Fuca Ridge, Bemis et al. [1993] estimated an average heat flux of 0.1-3.1 MW per vent. Ginster et al. [1994] indicated an average heat output of 3.1 MW for single high-temperature vents on the southern JDFR and higher values on the ridge's northerly Endeavor segment, where hydrothermal plumes can reach 300 m above the seafloor. There their estimated heat flux density of 0.139 MW m would, over a source area of 10 × 10 m, give a Q very much like the 13 MW heat flux value used in these calculations. Schultz et al. [1992] estimated a heat flux of 58 MW from a 4 × 5 m sulfide mound at Endeavor.
Calculations were made for a site latitude of 45°N. With the nonvertical component of rotation set to zero, the rotation vector (, , ) was (0, 0, 5.14 × 10 s); an example of the effect of nonzero [Garwood et al., 1985] on a hydrothermal plume is given by Lavelle [1995]. Mixing parameters were typically AHMIN = 10 m s and AVMIN = 10-4 m s. CS was set to 0.2, and PR, the ratio of turbulent viscosity to diffusivity, was given the value 1. Mixing length, ls, was made equal to (xyz)1/3 [e.g., Reynolds, 1990]. Values for lS, CS, and PR are not uniquely assignable, but values here are in ranges commonly used. Sensitivity of results to these parameters can be found, in part, in studies of axisymmetric plumes by Lavelle and Baker [1994] and of plumes from line segment sources by Lavelle [1995].
Time and length scales pertinent to the development of convective plumes from point and extended sources in rotating but otherwise quiescent background environments may be found in the recent work of Jones and Marshall [1993], Maxworthy and Narimousa [1994], and Speer and Marshall [1995]. For point sources and under stratified and depth-unlimited conditions, the external set of three variables, B, f, and N, allows two length scales to be derived. The first, well known from nonrotating conditions, is lN = (B/N)1/4 [e.g., Turner, 1973]. By replacing N by f, a second length scale is uncovered: lf = (B/f 3)1/4 [e.g., Speer and Marshall, 1995] or lROT = (B/f 3)½, where B in the last equation is buoyancy flux density rather than total buoyancy flux [Maxworthy and Narimousa, 1994]. Atmospheric and laboratory observations for plumes rising in nonrotating environments without cross flow have established that maximum rise height, hRISE, scales like lN: hRISE ~ 3.75 lN [e.g., Hanna et al., 1982]. An example in which lROT was found to be relevant comes from rotating tank studies of Maxworthy and Narimousa [1994] on brine convection. They found that once a brine from a distributed (i.e., nonpoint) source reached a fluid depth of the order of 12.7 lROT, convection transitioned from a three-dimensional turbulent condition to one dominated by descending vortical columns. Length scales such as l and lROT help organize experimental, observational, and occasionally numerical [e.g., Speer and Marshall, 1995; Lavelle and Smith, 1996] results, but in all cases the coefficient of scaling must be determined from data.
Previous modeling work on hydrothermal megaplumes, which are caused by short-lived thermal releases associated with episodic tectonic events [e.g., Embley et al., 1995], have confirmed the utility of those scales in such cases. Speer and Marshall [1995] used lf and the timescale for Coriolis effects to be important, f -1, to scale rotational velocities in the counterrotating vortices that should develop during megaplume formation. Lavelle and Baker [1994] and Speer and Marshall [1995] found that l scales rise height as it does in the nonrotating case under typical benthic ocean conditions, and Lavelle [1995] has shown that the timescale for megaplumes reaching hRISE is ~ 4N. The fact that megaplumes form quickly permits the modeling assumption of a quiescent background.
In the case of chronically discharging hydrothermal plumes, the assumption of a quiescent background environment is no longer tenable. Furthermore, these more typical hydrothermal plumes are not products of ephemeral sources of heat. They do not produce bottom-detached lenses of water, the geostrophic adjustment of which Gill [1981] and McWilliams [1988] have studied. Consequently, one must not expect that such plumes are associated with a Burger number ~1 nor expect that the description of bent-over plumes should involve length scales l and lf, as is made clear below.
In the case of continuous discharging plumes in a cross flow, the set of external variables from which length scales can be constructed is increased from three to four: the variable U must be added to the previous set. Generalized length scales that result from the expanded set are B(1+) N U (-3-4) and B(1+) f U(-3-4). When = -1, advective length scales LNADV = U/N and LfADV = U/f emerge. When = -3/4, lN and lf are recovered. But can take other values. Atmospheric and laboratory observations of plume in a cross flow show that lN (i.e., = -1) is not an appropriate scaling for plume rise height, for example. Those data (and entrainment theory, e.g., Middleton [1986]) point to a value for that is more nearly -2/3, so that rise heights of plumes in cross flows scale like lCROSS = [B/(UN)]1/3 [e.g, Hanna et al., 1982], not like [BN]1/4 as in cases without background flow. For plumes in a cross flow, it is only through observations that the coefficient of scaling and indeed even the value of can be found. Neither l or lf have demonstrated roles in the scaling of cross-flow plume results.
Transverse width of the plume is an important measure of the likelihood of counterrotating vortex pair development seen in studies of plume formation in quiescent backgrounds [e.g., Speer, 1989]. But it is clear that the width of a plume transverse to the cross-flow direction cannot scale like LNADV, LfADV, l or lf, except in the limit U 0: plumes having the same source B grow narrower with increasing U, but none of those four scales permit plumes with that U0 dependence. The width must be inversely proportional to some power of cross-flow velocity, perhaps scaling as lCROSS. The same arguments can be made that l and lf play no role in setting plume height. For the cross-flow magnitudes studied in this paper, plume width above the source is always much less than lf, the size at which one might expect the Coriolis force to cause counterrotating vortices in the manner previously predicted for plumes without background flow [Lavelle and Baker, 1994; Speer and Marshall, 1995]. Balancing vertical and horizontal mass flux at the top of the plume stem in simple calculations leads to that conclusion too. For the smallest U0 studied in the paper, cross flow overwhelms upstream density-driven flow in the plume cap above the source. With U larger than the density-driven upstream velocity in the plume cap, the plume bends forward. With U and B typical of chronic hydrothermal discharge conditions, the cross-flow carries plume mass injected to the level of neutral buoyancy away at a rate that prevents a cross-flow dimension of the plume to approach lf there.
An interesting transition to the axisymmetric plume case should occur in a sequence of plumes as U 0. Then one can expect the plume to first develop an anvil shape, with some upstream progression of the plume and some widening of the plume in the transverse direction. At even smaller U, the plume must begin to take on all the conditions of plumes described in earlier papers, including a counterrotating vortex pair. U, of course, drops out of the scaling arguments in this limit and l and lf prevail. The study of that transition from cross flow to axisymmetric plumes has not been undertaken, nor have enough cross-flow experiments been run that the process of scaling confirmation analogous to those of Speer and Marshall [1995] or Lavelle and Smith [1996] is possible now.
Momentum equations were centered differenced in space in energy conservation form and leapfrogged in time. An Asselin [Asselin, 1972] filter having = 0.15 was applied each time step to control temporal mode splitting. Integrations involved solving a 3-D Poisson equation for p [Harlow and Welch, 1965] on a staggered grid each time step using direct method solvers HS3CRI and HS3CRT (R. Sweet, National Bureau of Standards, Boulder, Colorado, 1985). Transport equations were forward time differenced. Their diffusion terms were centered differenced in space, while advection terms were upstream differenced but corrected each time step for unwanted numerical diffusion using the procedure of Smolarkiewicz [1983] and Smolarkiewicz and Clark [1986]. Even in the presence of strong advection, upstream differencing with correction can maintain relatively sharp property gradients like those found at plume stem walls, while maintaining positive definiteness of calculated quantities. The last attribute proves useful in eliminating the occurrence of unphysically low property anomalies.
The accuracy of the model has been examined in the following ways. Requisite conservation of mass, momentum, heat, salt, and energy was carefully checked under a variety of forcing, as was symmetry (or antisymmetry) of all field variables under symmetric forcing. Conditions for the independence of results from model grid size were examined by Lavelle and Baker [1994] using the axisymmetric realization of the model. Results on sensitivity to the Smagorinsky coefficient and Prandtl number and on the model relationship of plume rise height to buoyancy flux B0 are found in the same place. Additional sensitivity experiments, specifically for rise height dependence on buoyancy flux and buoyancy frequency N, were conducted for a line segment hydrothermal source rising into a quiescent background environment [Lavelle, 1995]. That analysis demonstrated that the rise height of the model scales as B1/3 in a quiescent background setting, as would be expected for that source configuration. Model results on the scaling of plume rise height with cross-flow strength are also encouraging. Atmospheric observations show plume rise height hRISE depending on the cross-flow velocity U as U-, where is ~1/3 [Hanna et al., 1982]. As demonstrated later, this model shows a rise height dependence on cross-flow velocity of U-0.4, with a coefficient of proportionality of similar magnitude to atmospheric cases.
Additional comparisons of model results to laboratory and field observations have been made. The utility of the Smagorinksy-Lilly subgrid-scale turbulence formulation was examined [Lavelle and Smith, 1996] by comparing model with laboratory results [Fernando and Ching, 1993] for a brine convection in a rotating tank. The importance of self-generated turbulence during convection has long been recognized [e.g., Priestly, 1956], and those numerical experiments in conjunction with laboratory results have reconfirmed the view that time- and space-dependent turbulence closure is essential. Furthermore, this model, under conditions of cross flow in a nonrotating environment, produces flow patterns comparable to the numerical results of Sykes et al. [1986] and to experimental results they cite, as will be described more fully below. Finally, hydrographic profiles from this model for regions affected by chronic hydrothermal sources have the form and magnitude of perturbed hydrographic profiles observed in the field (J. W. Lavelle et al., Effects of deep ocean hydrothermal discharge on near-source hydrography: Surrogate field studies with a convective plume model, Pacific Marine Environmental Laboratory contribution 1735, National Oceanic and Atmospheric Administration, Seattle, Washington, 1996).
Model content vis-a-vis the simpler entrainment model of Morton et al. [1956] was examined by computing the effective entrainment coefficient from axisymmetric convective plume results [Lavelle and Baker, 1994]. Horizontal velocities at the stem wall were divided by the average vertical velocity at the same height in the stem, and results were profiled as a function of height from the buoyancy source. Entrainment profiles so derived showed the entrainment coefficient having a magnitude comparable to that customarily used with integral theory ( ~ 0.1). Contrary to the assumption of entrainment theory the entrainment coefficient was not constant with height, even becoming negative in the plume cap region [Lavelle and Baker, 1994]. The number of model checks and comparisons to measurements enumerated here demonstrates how well tested this convection model is.
Return to previous section or go to next section