Upper Ocean Physics as Relevant to Ecosystem Dynamics:

Abstract

Our picture of the carbon cycle in the upper ocean relies on ecosystem dynamics for an understanding of the efficiency of nutrient uptake and export of carbon in the form sinking carbon particles, and for the fraction of recycled and exported particulate carbon production. A fundamental variable appears to be the size distribution of phytoplankton. Also, ultimately, an understanding of the partitioning between calcitic and siliceous based ecosystems may be important to predicting the long-term ocean carbon cycle.

Ecosystem dynamics of the upper ocean are driven by the interplay between dynamics of the surface ocean mixed layer and the depth of light penetration. This interplay is illustrated by noting the effect of high-frequency fluctuations in mixed layer depth from a physical model (Archer, et al., 1993) on an ecosystem model developed at weathership station Papa in the subarctic Pacific ocean (Frost, 1987). Three families of surface ocean mixed layer model are available for use by plankton ecologists, and although the physical mechanisms by which mixing occurs differ among the model groups, all are in general successful at predicting the observed ocean mixed layer depth. This paper explores behavioral distinctions between the three types of models, and summarizes previously published comparisons of the generality, accuracy, and computational requirements of the three models. Nutrients are supplied to the euphotic zone by the exchange of water between nutrient-depleted surface waters and nutrient-rich deeper waters. Current understanding of this process is still problematic, with rates of mixing required to balance nutrient uptake estimates higher than values predicted based on turbulence studies. Evidence is reviewed that episodic mixing, driven by frontal and mesocale motions, may be responsible for a significant fraction of vertical nutrient transport.

Key Words: Ocean carbon cycle, CO₂

1. The Role of Plankton Ecology in Climate Change Research

1.1. The Biological Pump and Seawater Carbonate Chemistry

1.1.1. A Summary of Carbonate Buffer System Chemistry

1.1.2. Aqueous Chemistry

where K is a thermodynamic equilibrium constant. [Note: the activity of a solute represents the apparent concentration, which approaches the actual concentration near the limit of infinite dilution. Activity is denoted by braces (e.g., {CO₂}), whereas concentrations are indicated by brackets (as [CO₂]). A thermodynamic equilibrium constant is defined in terms of species activities, as above.]

Analytical determination of the relative activities of H₂CO₃ and CO₂ is difficult, so it is convenient to group the two species together using the notation:

H₂CO₃^* = H₂CO₃ + CO₂. (2)

Carbonic acid undergoes de-protonation to form bicarbonate and carbonate ion. The equilibrium expression for the first pH reaction is expressed using the H₂CO₃^* notation (this is called a composite equilibrium expression):

The second dissociation reaction is described by

HCO₃^- = H⁺ + CO₃⁼,

K₂ = F({H⁺}{CO₃⁼},{HCO₃^-}). (4)

For convenience, apparent equilibrium constants, denoted by primes, can be defined in terms of concentrations of the carbonate species, rather than activities, as

Numerous authors have published values for K^'₁ and K^'₂, (Hansson, 1973, Mehrbach, et al., 1973, Dickson and Millero, 1987, Roy, et al., 1993) and the selection of the best pair is still a matter for research and debate (Millero, et al., 1993, Takahashi, et al., 1994). One set of expressions for K^'₁ and K^'₂ is given in the box.

The interaction between the carbonate species in seawater leads to the aforementioned counter-intuitive inverse behavior between CO₂ and CO₃⁼. When carbonate ion is removed from solution (such as by formation of solid CaCO₃), bicarbonate ion undergoes a disproportionation reaction (until the equilibrium expressions (21) are satisfied):

The net result of this transaction is to increase the pCO₂ of the water parcel as CO₃⁼ is removed. (Contrast this with the simple-minded expectation that removal of a carbon species (CO₃⁼) should reduce the concentrations of all other carbon species, the "falling tide lowers all boats" scenario). Natural variations in seawater pH cause the relative concentrations of CO₂ and CO₃⁼ to vary inversely with each other, like the seats on a teeter-totter, while the relative changes in the concentration of HCO₃^- are much smaller.

The carbonate buffer system as written above has two degrees of freedom; that is, two measurements are required to specify the system entirely. The "alkalinity" of a sample is defined as the amount of acid required to neutralize all of the weak bases in the solution (primarily bicarbonate and carbonate ion), and can be measured by pH titration (see (Stumm and Morgan, 1981)). Other practical measurements include the concentration of [[Sigma]]CO₂ [defined as CO₂ + HCO₃^- + CO₃⁼: measurement technique by (Johnson et al., 1985), and the pCO₂ (Takahashi, et al., 1982), and pH (Clayton and Byrne, 1993). Measuring pCO₂ directly is the most reliable way of determining air-sea pCO₂ differences in the field. However, for modeling CO₂ in the ocean, pCO₂ should be calculated from the alkalinity and [[Sigma]]CO₂, since these are conservative to mixing. A simple algorithm for calculating pCO₂ from alk and [[Sigma]]CO₂ is given in the box.

1.1.3. Gas Equilibrium Chemistry

where K^'_[[Eta]] is called the Henry's law constant. The inverse of the Henry's law constant is referred to as K^'₀. The pCO₂, or "partial pressure" of CO₂, is defined as the gas pressure exerted by CO₂ in the head space above a water sample at equilibrium. Seawater that is not in gas-exchange equilibrium with the atmosphere will have a pCO₂ that is different from the atmospheric value.

The Henry's law coefficient, K^'_H, varies strongly with temperature, tending to higher pCO₂ at higher temperature. The [CO₂]_aq is also affected by temperature dependence of K^'₁ and K^'₂ (for a full treatment of these effects, see (Broecker and Peng, 1982), page 150). The net temperature dependence of seawater pCO₂ is an increase by about 4.23% per degree Celsius (Takahashi, et al., 1993).

The pCO₂ of sea water is buffered by the carbonate system, relative to a non-pH active gas like O₂. This is because a change in pCO₂ must be accompanied by changes in the concentrations of bicarbonate and carbonate ions. At the normal pH of seawater, the carbonate buffer system provides an effective reservoir for CO₂ that is roughly 10 times larger than the reservoir for an unbuffered gas like O₂. Therefore, the equilibration time of seawater pCO₂ is roughly 10 times slower than that for O₂. The "buffering strength" of seawater is often expressed as the Revelle buffer factor ( ), which is defined as the ratio between the relative sensitivity of pCO₂ and total CO₂ to an incremental addition or removal of CO₂ to a water parcel.

For example, the average value of in global ocean surface waters is about 10. Therefore, a 1% increase in TCO₂, about 20 umol kg^-1 in typical surface seawater with TCO₂ = 2000 umol kg^-1, will result in an increase in pCO₂ of 10%, which corresponds to an increase in H₂CO₃^* of approximately 2 umol kg^-1. The numerical correspondence between the factor of 10 decrease in gas equilibration response time described in the last paragraph and the value of the Revelle buffer factor of approximately 10 is merely a confusing coincidence. Waters with a smaller value of the Revelle factor can take up more CO₂. Global surface ocean waters range from 8 in low CO₂ tropical surface waters to 12 in high CO₂ polar waters, and in deep water containing high CO₂ concentrations reach as high as 17 (Takahashi, et al., 1980, Takahashi, et al., 1993).

1.2. Nutrient Limitation and the High Nutrient Low Productivity (HNLP) Problem

In many regions of the world's oceans (the North Pacific, the Southern Ocean, and the equatorial Pacific) surface production does not fully deplete the surface nutrients during the summer growing season. In these regions, higher efficiency of nutrient utilization (i.e. lower nutrient concentrations) could drive sea surface pCO₂ to lower values. Simple box models of the ocean - atmosphere system have shown that the atmospheric pCO₂ is particularly sensitive to nutrient concentration in the high-latitude oceans (Sarmiento and Toggweiler, 1984, Knox and McElroy, 1984, Siegenthaler and Wenk, 1984),. If nutrients in the high-latitude surface waters were completely utilized by primary production, then the pre-industrial atmospheric pCO₂ could have been as low as 165 ppm; lower efficiency of nutrient use could have increased pCO₂ to as high as 425 ppm (pre-industrial atmospheric pCO₂ was 280 ppm). Obviously, understanding the factors limiting growth in these regions, and predicting how they might respond to changing climate, is crucial to global CO₂ modeling. Unfortunately, the limitation of phytoplankton growth is still a subject of research and debate.

One hypothesis is that the rate of primary production is limited by grazing by microzooplankton (Evans and Parslow, 1985, Frost, 1987, Fenchel, 1988). Microzooplankton have short reproductive cycles, and so can respond quickly to variations in primary production, in contrast to macrozooplankton, which have annual reproductive cycles. In the North Atlantic, where the spring bloom depletes the surface nutrients, it could be that the phytoplankton escape predation in the spring because of the deeper wintertime mixed layer (which depletes microzooplankton stocks to low levels (Frost, 1987)).

Another hypothesis involves the micronutrient iron, which may be limiting the primary production rate in some regions (Martin and Gordon, 1988, Martin and Fitzwater, 1988). Iron is thought to be used by phytoplankton in the enzymatic system for conversion of nitrate to ammonia; thus, a limitation of iron would inhibit new (nitrate-based) growth (Reuter and Ades, 1987, Morel, et al., 1991). Growth experiments using phytoplankton samples from the North Pacific were performed both with and without the addition of trace quantities of iron, and in most cases, the iron-added system showed enhanced phytoplankton growth rates relative to the low-Fe blank. Iron is supplied to the surface ocean by atmospheric flux of continental dust (Martin and Gordon, 1988), and regions of low nutrient utilization tend to correspond to regions of low atmospheric particulate input.

Recent work in the subarctic Pacific (the Subarctic Pacific Ecosystem Research program) has generated what has become known as the "SUPER synthesis" of these two hypotheses (Miller, et al., 1991). The growth rates of small phytoplankton appear to be unlimited by the availability of iron; that is, in in situ conditions, the doubling time of picophytoplankton approach physiological limits. However, production of picoplankton is limited by microbial grazing (Frost, 1987). When iron is added to these waters, after a few days nutrient uptake becomes dominated by large diatoms (Martin and Fitzwater, 1988), which are able to escape predation and bloom to the limit of available macronutrients. Thus, understanding the factors that control phytoplankton size appears to be of fundamental importance to a predictive understanding of the HNLP phenomenon. The efficiency of nutrient recycling must also be closely controlled by plankton size, with lower efficiency in systems dominated by larger grazers which produce quickly sinking fecal material, relative to the slow sinking velocity of microbial grazer fecal material. In addition, the potential for fish production is also fundamentally limited by the size characteristics of organisms at the bottom of the food chain (Ryther, 1969). Thus a predictive model for the size distribution of plankton as a function of environmental conditions would be extremely useful to prediction of air-sea carbon fluxes, nutrient utilization efficiency, new vs. recycled production, and fisheries yields.

1.3. Partitioning between Calcitic and Siliceous Ecosystems

In addition to their alleged long-term effect on atmospheric CO₂, the functional characteristics of phytoplankton production can also be significant to trace gases such as dimethyl sulfide, produced by coccolithophorids (Keller, et al., 1989). DMS production is one of the major terms in the natural sulfur cycle of the atmosphere, and has been implicated in the formation of cloud condensation nuclii (Charlson, et al., 1987).

2. The Role of Upper Ocean Physics in Plankton Ecosystem Modeling

2.1. The Depth of Turbulent Mixing and Light Penetration

We illustrate the impact of mixed layer depth fluctuations on plankton population dynamics using a high-resolution one-dimensional coupled physics / ecosystem model of the mixed layer from weathership station Papa in the subarctic Pacific. The physical simulation was described by (Archer, et al., 1993), and ground-truthed to observations from the Canadian Weathership Program during the years 1973-1978. Model temperature structure is compared with observations in Fig. 1. The biological ecosystem model was developed as part of the Subarctic Pacific Ecosystem Research program (Frost, 1987, Miller, et al., 1991). The ecosystem model was run off-line; i.e. daily values of mixed layer depths and sea surface temperatures from the model, and sunlight intensity data from the weathership, were used to drive the ecosystem model. The parameters in the ecosystem model are appropriate to conditions of relatively low mixed-layer depth variability (Frost, 1987), and the response of the model most resembles data from Weathership Station Papa in the low-variability simulations.

Springtime warming drives shoaling of the mixed layer, which confines phytoplankton to the well-lit depths near the surface. Thus, fluctuations in mixed layer depth drive variability in the chlorophyll concentration of the model (Figure 2). When the model mixed layer depths are sub-sampled on longer time intervals, the variability of the chlorophyll concentration decreases (Figure 3). This example illustrates how upper ocean mixing dynamics can drive the behavior of the planktonic ecosystem, and motivates the coupling of ecosystem dynamics models with models of physical mixing. The next section explains the inputs required by mixing models, and reviews the existing literature on the physics of mixing in the upper ocean.

2.2. Air-Sea Fluxes

2.2.1. Heat

The "bulk" formulas for sensible (H) and latent (L_v E) heat fluxes involve drag (transfer) coefficients (C_H for latent and C_E for evaporative), the wind speed at a height of 10 meters (U₁₀), and the difference in temperature (T, where the subscript "s" denotes water surface, and "a" denotes air) or water saturation (relative humidity, q), respectively (from (Gill, 1982)):

H = [[rho]] _airc_p C_H U₁₀ (T_s - T_a) (8)

and

L_v E = [[rho]]_air L_v C_E U₁₀ (q_s - q_a). (9)

The best value for C_H ranges from 0.83 * 10^-3 (stable atmospheric conditions) to 1.10 * 10^-3 (unstable conditions). C_E is estimated to be 1.5 * 10^-3.

The formulation for the heat flux caused by radiative transfer is somewhat more complicated. The input radiation (Q_I,0) can be described as a function of latitude and season (Frouin, et al., 1989). Some of the incoming radiation is blocked by clouds in the atmosphere. The effect of clouds represents a major uncertainty in the heat balance of the surface ocean (Seager, et al., 1988, Wigley, 1989, Mitchell, et al., 1989). In addition to oceanic heat gain by incoming light, some heat is lost by outgoing blackbody radiation. There exist several formulations to calculate the net radiative heat flux; the following was presented by (Gill, 1982):

Q_I = Q_I,0 (1 - [[alpha]]_s) (1 - 0.7n_c) (10)

and

Q_B = 0.985 [[sigma]] T_s⁴ (0.39 - 0.05 e_a^1/2) (1 - 0.6n_c²), (11)

where Q_I is the downward radiative heat flux, Q_b is the upward long-wave radiation, [[alpha]]_s is the surface albedo (roughly 0.3), n_c is the fractional cloud cover, [[sigma]] is the Boltzmann constant, T_s is the sea surface temperature, and e_a is the vapor pressure of the atmosphere (mbar). The standard deviation of the uncertainty of the total heat flux has been estimated to be on the order of 25 W m^-2 (Weare, et al., 1981); see also (Seager, et al., 1988, Davis, et al., 1981). For comparison, the global average heat flux is on the order of 240 W m^-2 (Gill, 1982), and the perturbation due to the anthropogenic CO₂ will be approximately 4 W m^-2 (Lorius, et al., 1991).

2.2.2. Fresh Water

2.2.3. Light

Q(z) = Q_blue,I e^{-k' z} (12)

and

k' = .04 + .0088 Chl + .054 Chl^2/3 (13)

where Q_blue,I is the incidence of blue visible light at the sea surface (roughly half of the total radiative energy), and Chl is the concentration of chlorophyll, and a and b are constants (from (Parsons, et al., 1984)). Only the blue, penetrating fraction of the light is photosynthetically active radiation (PAR). The rest of the light energy (Q_red,I) is red light, which is absorbed in the top meter. State-of-the-art models of PAR and light energy input to the mixed layer include such factors as wavelength dependent absorption, latitude, season, and angle of solar input (Sathyendranath and Platt, 1988, Morel, 1988).

2.2.4. Momentum

[[tau]] = [[rho]]_air C_D U₁₀². (14)

The drag coefficient, C_D is often represented by a constant. However, the high seas which result from heavy winds tend to increase C_D. Values of the drag coefficient as a function of wind speed were summarized by (Wu, 1982):

C₁₀ = 10^-3 (0.8 + 0.065 U₁₀) (15)

When this expression is incorporated into (7), the momentum flux into the surface layer scales as the wind speed to the third power. Since the momentum flux is a non-linear function of the wind speed, the use of wind speeds averaged over time scales of weeks to months results in underestimation of wind-driven mixing (for example, see (Rosati and Miyakoda, 1988))

2.3. Mixing Models

2.3.1. Integrated TKE, or "Bulk", Models

where

[[Lambda]] = 0 if f(dh,dt) < 0

and

T_s is the sea surface temperature, and [[Delta]]T is the temperature difference across the mixed layer base.

The two degrees of freedom in equation 16 are the sea surface temperature and the depth of the mixed layer; another independent equation is necessary to close the system. The requisite constraint is based on a kinetic energy budget. The source of TKE is the wind stress (a scaling coefficient, m, times the "friction velocity" of the wind, U_*³) and the sink is dissipation (D). The residual kinetic energy (the difference between input and dissipation) is transformed into potential energy by mechanical mixing of cold water into the mixed layer (also, the average TKE of the mixed layer decreases as quieter water is entrained):

Since (Kraus and Turner, 1967), many other workers have used models based on this formulation, and much of the work has centered on the dissipation of turbulence (see (Stevenson, 1979) for a review). The generation and dissipation of turbulence is the greatest weakness of this type of model; in the case of dissipation in particular, the link between the model and the actual physics is weak. Some authors scale TKE dissipation as a constant loss rate times the mixed layer depth (Kim, 1976, Niiler and Kraus, 1977). Others calculate the dissipation rate as a function of the total amount of mixed-layer TKE, the coriolis parameter (Garwood, 1977) and/or the Monin-Obukhov length scale for turbulence (Gaspar, 1988). To simulate the mixed layer response to hurricanes, (Elsberry, et al., 1976) used a parameterization for dissipation as a fraction of surface input of TKE that increases with increasing mixed layer depth.

The other weakness of the bulk TKE family of models is the generation of TKE (calculated as a constant, m, times the wind stress). The value of m is "tuned" to fit the model predictions to observed data, and the best value for m tends to vary with location and conditions (ranging from 0.1 to 0.39 (Martin, 1985) 0.3 to 0.9 (Price, et al., 1978) 0.4 to 0.5 (Davis, et al., 1981)). The instability of m and the uncertainty about TKE dissipation certainly detract from the predictive ability of the bulk TKE-type model.

2.3.2. Shear Instability Models

2.3.2.1. Formulation

and is called the gradient Richardson number.

The Richardson number used in most mixed layer models is defined somewhat differently. The mixed layer is viewed as a slab, with uniform velocity and density. Instead of the differential quantities [[partialdiff]][[rho]]/[[partialdiff]]z and [[partialdiff]]u/[[partialdiff]]z, the differences [[Delta]][[rho]]/h and [[Delta]]u/h are used, where [[Delta]][[rho]] and [[Delta]]u are the changes in density and velocity across the mixed layer base. The length scale is the thickness of the mixed layer, h, rather than a scale associated with the thickness of the mixed layer base. The dimensionless number thus defined,

is called the bulk Richardson number (Pollard, et al., 1973, Price, et al., 1986). In the shear models, mixing is assumed to begin when Ri_b falls below a critical value, and water is entrained until Ri_b reaches the critical value once again. Models based on the bulk Richardson number predict mixed layer variations quite well, and the only potentially "tunable" parameter is the critical Richardson number, which appears to be stable throughout a variety of locations and climatic conditions (Price, et al., 1978, Price, et al., 1986). Theoretical justification for use of the bulk, rather than the gradient, Richardson number for these models is still a matter of discussion (Pollard, et al., 1973, Price, et al., 1978).

The shear instability model of (Price, et al., 1986) uses both Ri_b and Ri_g to determine mixing (Figure 6). The model was used to simulate detailed diurnal fluctuations in the mixed layer depth and current profiles from the subtropical Pacific. In the model, the heat fluxes (except penetrating solar radiation) and the wind stress are applied to the mixed layer. The bulk Richardson number at the base of the surface box is calculated, and if the shear is greater than the density stratification necessary to support it (Ri_b < 0.65), then the properties of the two adjacent boxes are averaged (the boxes are mixed). This process continues downward until Ri_b > 0.65. The mixed layer is considered homogeneous for momentum and density, as it is in the bulk TKE formulations considered above. The density transition at the mixed layer base is smoothed using the gradient Richardson number. If Ri_g between two adjacent boxes is smaller than a critical value (Ri_g < 0.25), the adjacent boxes are partially mixed until the shear becomes sub-critical again, and the process is repeated until Ri_g is greater than or equal to the critical value throughout the water column. Thus, there is a region of partial mixing below the completely mixed zone. The model also includes convective entrainment driven by surface cooling.

2.3.2.2. Behavior

Also, the value of the critical Ri_b required to simulate oceanographic data using a shear instability mixed layer model is fairly constant throughout a range of climatic conditions (Price, et al., 1978, Pollard, et al., 1973, Price, et al., 1986). This can be contrasted with the wind-scaling coefficient used in the TKE model, which must be tuned to simulate data from different locations or climatic regimes.

2.3.3. Turbulence Closure Models

2.3.3.1. Formulation

and

where U is velocity, [[Theta]] is temperature, f is the coriolis parameter, P is pressure, g is the gravitational constant, [[nu]] is the molecular viscosity, [[alpha]] is the thermal diffusivity, and [[beta]] is the thermal expansion coefficient. The equations for the mean flow and temperature can be found by decomposing the total velocity and temperature fields into mean and fluctuating components and averaging:

and

The terms xto(u u) and xto(u [[theta]]) represent the effect of the turbulent fluctuations on the mean velocity and temperature (the Reynold's stress and "eddy diffusion"). Finding analytical expressions for these terms is the problem at the heart of the turbulence closure models.

Derivation of expressions for these terms begins with these two sets of equations, and incorporates several empirical simplifying relationships based on laboratory data (the interested reader is referred to (Mellor, 1973)). The salient point to understand about the derivation is that the solutions contain several empirical constants, which are fit to laboratory data from neutral (unstratified) flows. (Mellor and Yamada, 1982) show the surprising result that the model, using these constants, is also able to predict laboratory data for stratified flow, flow in a pipe, and flow over a curved surface. Some of the terms in the equations also require a constant with units of length to maintain dimensional homogeneity. Mellor assumes that the length scales are all proportional to a single "mixing length scale" (l). This assumption, and the derivation of a number for l, is considered a weak link in the turbulence closure scheme. The choice of l is crucial; the levels 2 and 2.5 models (described below) supply eddy diffusion coefficients that scale directly to l.

The full version of the model is somewhat cumbersome, and some of the terms can be neglected without loss of predictive value. (Mellor and Yamada, 1974) presents a hierarchy of closure schemes which employ increasing levels of simplification, based on the dependence of each term on the degree of anisotrophy (directional inhomogeneity) of the turbulence (the term is called "a", and is typically around 15% for the ocean surface boundary layer). The level 4 model is simplified the least, and is in general too complicated for routine oceanographic use, without any corresponding benefit in predictive value. The level 3 model neglects all terms which clearly scale with a². There is some uncertainty to this scaling procedure. For the level 3 model the advective and diffusive terms for q² (a synonym for 2 * TKE) are assumed to scale with a¹ and are left in; for the level 2 model, they are assumed to scale with a² and are neglected. This has the importance consequence in the level 2 model of setting the diffusive flux of TKE equal to zero, so that all turbulence is dissipated locally to its generation (Niiler and Kraus, 1977).

Mellor's model of choice is called the level 2.5 model (Mellor and Yamada, 1982). Level 2.5 starts from the level 3 model but neglects the advection and diffusion of variance of the material terms like temperature and salinity (equations for
xto([[theta]]²) and xto(s²)). The result is the elimination of several differential equations; the savings depends on the number of scalar quantities being modeled (e.g., nutrients, pCO₂).

The final, reduced form of the level 2.5 model adds two equations to the averaged equations of motion that would be required to model any geophysical fluid. The first describes the time and space distribution of q²/2, the turbulent kinetic energy. Energy is added at the top as a boundary condition, advects, and diffuses according to the calculated eddy diffusion coefficient. Turbulent energy is generated by shear at the expense of mean flow. The length scale, l, is determined using a second differential equation for the quantity q²l. The result of this formulation is that l, the length scale of the largest eddies, is time dependent, has memory, and is allowed to be transported. The eddy diffusion coefficient at a given time and depth is calculated as a function of q² and l. Once the diffusion coefficient field is known, it is used for mixing the scalar quantities.

The mixing coefficient predicted by the model turns out to be a function of the local gradient Richardson number (Niiler and Kraus, 1977, Mellor and Durbin, 1975); above Ri_g = 0.23, no mixing at all is predicted. This is similar to the Price model, in which mixing takes place until a critical Ri_g > 0.25 is attained. The level 2.5 model has been used in a global ocean model that focuses on the surface ocean and the sea surface temperature (Rosati and Miyakoda, 1988).

2.3.3.2. Behavior

One obvious distinction is the degree of mixing in the turbulence closure model mixed layer. Whereas the other formulations assume homogeneity of all properties within the mixed layer, the turbulence models predict high but finite mixing coefficients within this zone. Under conditions of deep convective mixing, the "bulk" assumption (that the mixed layer is completely homogeneous in all properties) is probably inadequate, but under conditions of shallower, wind-driven mixing, the "bulk" approximation may be closer to reality than the output of the turbulence models (Martin, 1986).

The fundamental distinction between bulk and shear models presented earlier was the source of the turbulent energy for entrainment: from wind stress directly or from current shear. The main source of entrainment energy in the turbulence closure models is difficult to judge from first principles. For the primary source of turbulence to be generation at the surface, the rate of transport of turbulent energy through the mixed layer must be greater than the rate of generation of turbulence by shear at the mixed layer base. In the level 2 model, transport of turbulence does not occur; all turbulence is dissipated locally. Numerical experiments using levels 2.5 and higher have shown that the transport flux of turbulent energy from the surface zone in these models is also small relative to the generation of turbulence by shear (Klein and Coantic, 1981, Martin, 1986). Thus the turbulence closure models appear to function primarily as shear instability models.

2.3.4. Mixed Layer Model Comparison

2.3.4.1. The Stability of Empirical Parameters

2.3.4.2. Accuracy

2.3.4.3. Computational Efficiency

2.4. Physical Transport of Dissolved Nutrients into the Euphotic Zone

A simple model of nutrient cycling in the upper ocean was presented by (Peng, et al., 1987), to reproduce data from the Greenland sill. The upper box of the model is the mixed layer; the lower box represents the thermocline (Figure 7). The boundary between the boxes fluctuates according to the observed annual cycle of mixed layer depth. Both boxes are well-mixed. Organic carbon is produced in the upper box and regenerated within the lower box. When the mixed layer deepens in the winter, nutrients from the deep box are entrained into the mixed layer. The model neglects the diffusive flux of nutrients across the boundary (either between the two boxes or from the deep ocean into the lower box) and the sinking flux of particles at the bottom boundary. However, the annual cycle of surface water nutrient concentration and the dissolved O₂ and CO₂ predicted by this simple model are similar to the observations.

When the vertical density structure of the water column and the processes of nutrient transport and regeneration are considered in greater detail, the picture becomes more complex (Figure 8). Although most of the organic matter that falls from the euphotic zone is regenerated in the upper few hundreds of meters of the water column, some fraction (~25% at 200 meters (Martin, et al., 1987)) inevitably escapes from the deepest reach of the mixed layer; it is this fraction that drives respiration in the deep water and benthos. This loss must be balanced by upward diffusion or advection of nutrients.

In the subpolar gyres, the net upwelling driven by the wind stress curl is an obvious mechanism for transporting nutrients to the euphotic zone. In the subtropical gyres, the net vertical velocity is downward, and accounting for the nutrient supply inferred from the sinking carbon flux is more difficult. (Hayward, 1987) documented the problem in detail in the central North Pacific, using data from several California / Hawaii transects and quasi-time series data from a site near Hawaii. The surface nutrient structure in this region is characterized by low nutrient concentrations in the euphotic zone, on top of a "nutricline", or region of nutrient increase with depth. The top of the nutricline is located well below the depth of the seasonal thermocline, and is roughly at the depth of the deepest reach of the mixed layer in winter. No seasonal structure in the nutrient/depth relation was detected in the nutricline data, which would be indicative of entrainment of the nutrient rich thermocline water into the winter mixed layer. Calculation of nutrient mixing rates through the nutricline, based on the paradigm of Fickian eddy mixing, yielded a nutrient supply rate that was an order of magnitude too small to account for the observed production. (This calculation is sensitive to the appropriate value of the eddy mixing coefficient, however, which is not well known.) Perhaps more importantly, Hayward observes that there was no correlation between the nutrient gradient and the productivity rate, as would be predicted by a Fickian eddy mixing model. He concludes that nutrient supply to the euphotic zone must be more complicated than can be explained by the eddy mixing formulation. The nutrient supply problem has also been documented in the subtropical Atlantic near Bermuda (Jenkins, 1987), with the conclusion that the Fickian "background" eddy mixing rate is too slow by an order of magnitude to mix nutrients into the euphotic zone fast enough to meet the observed production demand.

Several potential solutions to the nutrient supply problem have been suggested. (Villareal, et al., 1993) proposed that diatoms in this region may be able to migrate vertically to reach both the deep nutrients and the shallow euphotic zone in the course of a 24-hour cycle). There are also abundant indications that nutrient transport may be an episodic process, rather than a "constant background" eddy diffusive flux. (Jenkins, 1987) observed an anomalous nutrient profile in the subtropical Atlantic that appeared to represent an injection of high nutrient thermocline water into the euphotic zone. The apparent transport of nitrate to the euphotic zone during this event was high enough to supply 20-30% of the annual production demand. Based on a calculated residence time of inorganic nitrate, Jenkins calculates that these episodic events could supply the entire nutrient requirement of the euphotic zone, but the probability of actually observing such an event on any given cast would only be around 1%. Abundant theoretical evidence exists that nutrient transport may be dominated by motions on a frontal and mesoscale in the upper ocean. (Bleck, et al., 1988) constructed a two dimensional isopycnal vertical coordinate numerical model of a convergent front. Using typical open ocean conditions, they calculated local time-averaged vertical velocities associated with fronts of greater than 5 meters day^-1. (Pollard and Regier, 1990) deployed a towed undulating CTD and an acoustic doppler current profiler in the subtropical Atlantic and inferred that the mechanisms predicted by (Bleck, et al., 1988) were correct, and that there were vertical velocities of tens of meters per day at that location. Intense fronts are found at the boundaries between water masses (gyres) in the upper ocean; less intense fronts also form during the decay of mesoscale eddies, which are thought to be ubiquitous in the ocean ((Semtner and Chervin, 1992) and references therein). (Strass, 1992) presented towed undulating vehicle data that included a fluorescence sensor for chlorophyll in addition to CTD measurements. He found the highest values of the subsurface chlorophyll concentration on the anticyclonic side of a front, just as predicted. Also, he found that the spatial scale of chlorophyll variability was similar to the Rossby radius of deformation, as predicted by (Woods, 1988). The similarity between the spatial scale of chlorophyll patchiness and the physical scale of eddy dynamics (the mesoscale) lends support to the idea that meso and frontal scale dynamics control the transport of nutrients to the euphotic zone in the subtropical gyre.

We conclude that while macro-nutrient supply (NO₃ and PO₄) clearly regulates photosynthesis in some pelagic systems, most notably in the subtropical gyres, the mechanisms by which nutrients are supplied to the euphotic zone are still a topic of active research and that an off-the-shelf nutrient transport model is as yet not available. In locations such as the subtropical Pacific, where nutrient supply clearly has a dominant role in determining phytoplankton growth and depth distribution (Hayward, 1987), a simple eddy diffusion formulation must be used, with the understanding that essential characteristics of nutrient transport, and thus their effect on phytoplankton dynamics, may be missing.

Summary

Citations

Archer, D. E. and E. Maier-Reimer. 1994. Effect of deep-sea sedimentary calcite preservation on atmospheric CO₂ concentration. Nature 367: 260-264.

Baumgartner, A. and E. Reichel. 1975. The world water balance: Mean annual global, continental and maritime precipitation, evaporation and run-off. Elsevier, Amsterdam, Holland.

Bleck, R., R. Onken, and J. D. Woods. 1988. A two-dimensional model of mesoscale frontogenesis in the ocean. Quarterly Journal of the Royal Meterological Society 114: 347-371.

Broecker, W. S. 1991. Keeping global change honest. Global Biogeochemical Cycles 5: 191-192.

Broecker, W. S. and T. H. Peng. 1982. Tracers in the sea. Eldigio Press, Palisades, NY.

Charlson, R., J. E. Lovelock, M. O. Andreae, and S. G. Warren. 1987. Oceanic phytoplankton, atmospheric sulfur, cloud albedo and climate. Nature 326: 655-661.

Clayton, T. D. and R. H. Byrne. 1993. Calibration of m-cresol purple on the total hydrogen ion concentration scale and its application to CO₂-system characteristics. Deep-Sea Research 40: 2115-2129.

Davis, R. E., R. deSzoeke, and P. Niiler. 1981. Variability in the upper ocean during MILE. Part I: The heat and momentum balances. Deep Sea Research 28: 1427-1451.

Davis, R. E., R. dsSzoeke, and P. Niiler. 1981. Variability in the upper ocean during MILE. Part II: Modeling the mixed layer response. Deep Sea Research 28: 1453-1475.

Dickson, A. G., and F. J. Millero. 1987. A comparison of the equilibrium constants for the dissociation of carbonic acid in seawater media. Deep-Sea Research 34: 1733-1743.

Ellison, T. H., and J. S. Turner. 1959. Turbulent entrainment in stratified flows. Journal of Fluid Mechanics 6: 423-448.

Elsberry, R. L., T. S. Friam, and R. N. J. Trapnell. 1976. A mixed layer model of the oceanic thermal response to hurricanes. Journal Geophysical Research 81: 1153-1162.

Evans, G. T., and J. S. Parslow. 1985. A model of annual plankton cycles. Biological Oceanography 3: 327-347.

Fenchel, T. 1988. Marine plankton food chains. Annual Review of Ecological Systems 19: 19-38.

Frost, B. W. 1987. Grazing control of phytoplankton stock in the open subarctic Pacific Ocean : a model assessing the role of mesozooplankton, particularly the large calanoid copepod Neocalanus spp. Marine Ecology Progress Series 39: 49-68.

Frouin, R., D. W. Linger, C. Gautier, K. S. Baker, and R. C. Smith. 1989. A simple analytical formula to compute clear sky total and photosynthetically available solar irradiance at the ocean surface. Journal Geophysical Research 94 : 9731-9743.

Gargett, A. E., T. B. Sanford, and T. R. Osborn. 1979. Surface mixing layers in the Sargasso Sea. Journal Physical Oceanography 9: 1090-1111.

Garwood, R. W. 1977. An ocean mixed layer model capable of simulating cyclic states. Journal Physical Oceanography 7: 455-468.

Gaspar, P. 1988. Modeling the seasonal cycle of the upper ocean. Journal Physical Oceanography 18: 161-180.

Gill, A. E. 1982. Atmosphere-Ocean Dynamics. Acedemic Press, Orlando, Florida.

Hansson, I. 1973. A new set of acidity constants for carbonic acid and boric acid in seawater. Deep-sea Research 20: 461-478.

Hayward, T. L. 1987. The nutrient distribution and primary production in the central North Pacific. Deep-sea Research 34: 1593-1627.

Holm-Hansen, O, E.W. Helbling, and D. Lubin. 1993. Ultraviolet radiation in Antarctica: inhibtion of primary production. Photochemistry and Photobiology 58: 567-571.

Jenkins, W. J. 1987. The use of anthropogenic tritium and helium-3 to study subtropical gyre ventilation and circulation. Philosophical Transactions of the Royal Soceity of London 325: 43-61.

Jerlov, N. G. 1976. Marine Optics. Elsevier, Amsterdam, Netherlands

Johansson, O., and M. Wedborg. 1982. On the evaluation of potentiometric titrations of seawater with hydrochloric acid. Oceanologica Acta 5: 209-218.

Johnson, K.M., A.E. King, and J. Sieburth. 1985. Coulometric [[Sigma]]CO₂ analyses for marine studies: an introduction. Marine Chemistry 16: 61-82.

Keller, M. D., W. K. Bellows, and R. R. L. Guillard. 1989. A survey of DMS production in twelve classes of marine phytoplankton, pages 167-182in Saltzman, E. and W. Cooper, editors. Biogenic sulfur in the Environment. American Chemical Society, Washington, D.C.

Kim, J.-W. 1976. A generalized bulk model of the oceanic mixed layer. Journal of Physical Oceanography 6: 686-695

Klein, P. and M. Coantic. 1981. A numerical study of turbulent processes in the marine upper layers. Journal of Physical Oceanography 11: 849-863.

Knox, F. and M. McElroy. 1984. Change in atmospheric CO₂: Influence of the marine biota at high latitude. Journal of Geophysical Research 89: 4629-4637.

Kraus, E. B. and J. S. Turner. 1967. A one-dimensional model of the seasonal thermocline II. The general theory and its consequences. Tellus 19: 98-105.

Lisitzin, A. P. 1967. Basic relationships in distribution of modern siliceous sediments and their connection with climatic zonation. International Geology Review 9: 631-652, 842-863, 1114-1130.

Lorius, C., J. Jouzel, D. Raynaud, J. Hansen, and H. L. Treut. 1991. The ice-core record: climate sensitivity and future greenhouse warming. Nature 347: 139-145.

Martin, J. H. and R. M. Gordon. 1988. Northeast Pacific iron distributions in relation to phytoplankton productivity. Deep-sea Research 35: 177-196.

Martin, J. H., G. A. Knauer, D. M. Karl, and W. M. Broenkow. 1987. VERTEX: carbon cycling in the northeast Pacific. Deep Sea Research 34: 267-285.

Martin, P. J. 1985. Simulation of the mixed layer at OWS November and Papa with several models. Journal of Geophysical Research 90: 903-916.

Martin, P. J. 1986. Testing and comparison of several mixed-layer models. Naval Ocean Research and Development Activity. Technical Report Number 143.

Mehrbach, C., C. Culberson, J. E. Hawley, and R. M. Pytkowicz. 1973. Measurement of the apparent dissociation constants of carbonic acid in seawater at atmospheric pressure. Limnology and Oceanography 18: 897-907.

Mellor, G. L. 1973. Analytic prediction of the properties of stratified planetary surface layers. Journal of the Atmospheric Sciences 30: 1061-1069.

Mellor, G. L. and P. A. Durbin. 1975. The structure and dynamics of the ocean surface mixed layer. Journal of Physical Oceanography 5: 718-728.

Mellor, G. L. and T. Yamada. 1974. A hierarchy of turbulence closure models for planetary boundary layers. Journal of the Atmospheric Sciences 31: 1791-1806.

Mellor, G. L. and T. Yamada. 1982. Development of a turbulence clusure model for geophysical fluid problems. Reviews of Geophysics and Space Physics 20: 851-875.

Miller, C. B., B. W. Frost, P. A. Wheeler, M. R. Landry, N. Welschmeyer, and T. M. Powell. 1991. Ecological dynamics in the subarctic Pacific, a possibly iron-limited ecosystem. Limnology and Oceanography 36: 1600-1615.

Millero, F. J., R. H. Byrne, R. Wanninkhof, R. Feely, T. Clayton, P. Murphy and M. Lamb. 1993. The internal consistency of CO₂ measurements in the equatorial Pacific. Marine Chemistry 44: 269-280.

Mitchell, J. F. B., C. A. Senior, and J. W. Ingram. 1989. CO₂ and climate: a missing feedback? Nature, 341: 132-134.

Morel, A. 1988. Optical modeling of the upper ocean in relation to its biogenous matter content (Case I waters). Journal of Geophysical Research, 93 : 10749-68.

Morel, F. M. M., J. G. Reuter, and N. M. Price. 1991. Iron nutrition of phytoplankton. Oceanography 4: 71-78.

Niiler, P. P. and E. B. Kraus. 1977. One-dimensional models of the upper ocean, pages 143-172 in Kraus, E.B. editor, Modelling and Prediction of the Upper Layers of the Ocean. Pergammon Press, Oxford, U.K.

Oey, L.-Y., G. L. Mellor and R. I. Hires. 1985. A three-dimensional simulation of the Hudson-Raritan estuary. Part I: description of the model and model simulations. Journal of Physical Oceanography 15: 1676-1692.

Owens, N. J. P., J. N. Galloway, and R. A. Duce. 1992. Episodic atmospheric nitrogen deposition to oligotrophic oceans. Nature 357: 397-399.

Parsons, T. R., M. Takahashi, and B. Hargrave. 1984. Biological Oceanographic Processes. Pergammon Press, Oxford, U.K.

Peixoto, J. P. and A. H. Oort. 1983. The atmospheric branch of the hydrological cycle and climate, in al., A. S.-P. e. Variations of the global water budget. Reidel, London, U.K.

Peixoto, J. P. and A. H. Oort. 1992. Physics of Climate. 520 Pages. American Institute of Physics, New York.

Peng, T.-H., T. Takahashi, W. S. Broecker, and J. Olafsson. 1987. Seasonal variability of carbon dioxide, nutrients, and oxygen in the northern North Atlantic surface water: Observations and a model. Tellus 39B: 439-458.

Pollard, R. T. and L. Regier. 1990. Large variations in potential vorticity at small spatial scales in the upper ocean. Nature 348: 227-229.

Pollard, R. T., P. B. Rhines, and R. O. R. Y. Thompson. 1973. The deepening of the wind-mixed layer. Geophysical Fluid Dynamics 3: 381-404.

Price, J. F., C. N. K. Mooers, and L. Van. 1978. Observations and simulation of storm-induced mixed-layer deepening. Journal of Physical Oceanography 8: 582-599.

Price, J. F., R. A. Weller. and R. Pinkle. 1986. Diurnal cycling: observations and models of the upper ocean response to diurnal heating, cooling and wind mixing. Journal of Geophysical Research 91: 8411-8427.

Redfield, A. C. 1934. On the proportions of organic derivatives in sea water and their relation to the composition of plankton, in James Johnson Memroial Vol. Liverpool, U.K.

Reuter, J. G. and D. R. Ades. 1987. The role of iron nutrition in photosynthesis and nitrogen assimilation in Scendesmus Quadricuada (Chlorophyceae). Journal of Phycology 23: 452-457.

Riebesell, U., D. A. Wolf-Gladrow, and V. Smetacek. 1993. Carbon dioxide limitation of marine phytoplankton growth rates. Nature 361: 249-251.

Rosati, A. and K. Miyakoda. 1988. A general circulation model for upper ocean simulation. Journal of Physical Oceanography 18: 1601-1626.

Roy, R. N., L. N. Roy, K. M. Vogel, C. Porter-Moore, T. Pearson, C. E. Goode, W. Davis, F. J. Millero, and D. M. Campbell. 1993. The dissociation constants of carbonic acid in seawater at salinities 5 to 45 and temperatures 0 to 45^o C. Marine Chemistry 44: 249-267.

Ryther, J. H. 1969. Photosynthesis and fish production in the sea. Science 166: 72-76.

Sarmiento, J. L. and R. Toggweiler. 1984. A new model for the role of the oceans in determining atmospheric pCO₂. Nature 308: 621-624.

Sathyendranath, S. and T. Platt. 1988. The spectral irradiance field at the surface and in the interior of the ocean: a model for applications in oceanography and remote sensing. Journal of Geophysical Research 93: 9207-9280.

Seager, R., S. E. Zebiak, and M. A. Cane. 1988. A model of the tropical pacific sea surface temperature climatoloty. Journal of Geophysical Research, 93: 1265-1280.

Semtner, A. J. and R. M. Chervin. 1992. Ocean general circulation from a global eddy-resolving model. Journal of Geophysical Research 97: 5493-5551.

Siegenthaler, U. and T. Wenk. 1984. Rapid atmospheric CO₂ variations and ocean circulation. Nature 308: 624-626

Smith, R.C., B.B. Prezelin, K.S. Baker, R.R. Bidigare, N.P. Boucher, T. Coley, D. Karentz, S. MacIntyre, H.A. Matlich, D. Menzies, M. Ondrusek, and Z. Wan. 1992. Ozone depletion: Ultraviolet radiation and phytoplankton biology in Antarctic waters. Science 255: 952-959.

Stevenson, J. W. 1979. On the effect of dissipation on seasonal thermocline models. Journal of Physical Oceanography 9: 57-64.

Strass, V. H. 1992. Chlorophyll patchiness caused by mesoscale upwelling at fronts. Deep-sea Research 39: 75-96.

Stumm, W. and J. J. Morgan. 1981. Aquatic Chemistry. Wiley, New York.

Sverdrup, H. U. 1953. On conditions for the vernal blooming of phytoplankton. J. Cons. Explor. Mer. 18: 287-295.

Takahashi, T., W. S. Broecker, A. E. Bainbridge, and R. F. Weiss. 1980. Carbonate chemistry of the Atlantic, Pacific, and Indian Oceans: The results of the GEOSECS expeditions, 1972-1978. National Science Foundation, Washington D.C.

Takahashi, T., W. S. Broecker ,and S. Langer. 1985. Redfield ratio based on chemical data from isopycnal surfaces. J. Geophys. Research 90: 6907-6924.

Takahashi, T., D. Chipman, N. Schechtman, J. Goddard, and R. Wannikof. 1982. Measurements of the partial pressure of CO₂ in discreet water samples during the North Atlantic Expedition, the Transient Tracers of Oceans Project. Technical Report to NSF. Lamont-Doherty Geological Observatory, Palisades, NY.

Takahashi, T., J. Olafsson, J. G. Goddard, D. W. Chipman, and S. C. Sutherland. 1993. Seasonal variation of CO₂ and Nutrients in the high-latitude surface oceans: a comparative study. Global Biogeochemical Cycles 7: 843-878.

Therry, G. and P. Lacarrere. 1983. Improving the eddy kinetic energy model for planteary boundary layer description. Boundary-Layer Meteorology 25: 63-88.

Villareal, T. A., M. A. Altabet, and K. Culver-Rymsza. 1993. Nitrogen transport by vertically migrating diatom mats in the North Pacific Ocean. Nature 363: 709-712.

Weare, B. C., P. T. Strub, and M. D. Samuel. 1981. Annual mean surface heat fluxes in the tropical Pacific Ocean. Journal of Physical Oceanography 11: 705-717.

Weiss, R. 1974. Carbon dioxide in water and seawater. The solubility of a non-ideal gas. Marine Chemistry 2: 203-215.

Wigley, T. M. L. 1989. Possible climate change due to SO2-derived cloud condensation nuclei. Nature 339: 365-367.

Woods, J. 1988. Mesoscale upwelling and primary production, pages 7-38 in Rothschild, B. J., editor, Toward a theory on biological-physical interactions in the world ocean. Kluwer, Amsterdam, Netherlands.

Woods, J. D. and W. Barkmann. 1986. The response of the upper ocean to solar heating. I: The mixed layer. Quarterly Journal of the Royal Meteorological Society 112: 1-42.

Wu, J. 1982. Wind-stress coeffecients over sea surface from breeze to hurricane. Journal of Geophysical Research 87: 9704-9706.

Box

The equilibrium constant for the process

CO₂(g) = H₂CO₃^*(aq)

is given by (Weiss, 1974) as

ln K^'₀ = 93.4517 (100/T) - 60.2409 + 23.3585 ln (T/100) +
S { 0.0023517 - 0.023656 (T/100) + 0.0047036 (T/100)²}

For the dissociation of carbonic acid, (Dickson and Millero, 1987) give

K^'₁ = 10 ** {-(6320.18/T - 126.3405 + 19.568 ln T +
( 19.894 - 840.39/T - 3.0189 ln T) S^(0.5) + 0.00679S ) }

and

K^'₂ = 10 ** { -(5143.69/T - 90.1833 + 14.613 ln T +
( 17.176 - 690.59/T - 2.6719 ln T) S^(0.5) + .02169 S ) }

For boric acid, (Johansson and Wedborg, 1982) give

K^'_B = 10**{ -(1030.5/T + 5.5076 - 0.015469S
+ 1.5339^.10^-4 S²) }

where the total borate concentration of seawater can be estimated to be

[[Sigma]]B = 4.106 ^.10^-4 S / 35.

(B) Method of calculation of carbonate system parameters given values of alkalinity, [[Sigma]]CO₂, and salinity, from (Takahashi, et al., 1980).

TALK = total alkalinity, 10^-6 eq kg^-1

[[Sigma]]CO₂ and [[Sigma]]B in 10^-6 mol kg^-1

C1 = K^'₁/2.0

C2 = 1.0 - 4.0*K^'₂/K^'₁

C4 = [[Sigma]]B * K^'_B

AHT = 0.7E-8

DO ICNT = 1, 100

A = TALK - C4/(K^'_B + AHT)

X = A / TCO2

AH1 = C1 / X * (1.0 - X + SQRT( 1.0 + C2 * X * (-2.+X)))

IF(0.5E-4 .GE. ABS(1. - AHT-AH1)) EXIT

AHT = AH1

ENDDO

pCO2 = A * (AH1/K^'₁) / (K^'₀ * (1.0+2.0*K^'₂/AH1))