D. Archer1, S. Emerson2, T. Powell3, and C.S. Wong4
1Lamont-Doherty Earth Observatory of Columbia University
Palisades, NY
2School of Oceanography
University of Washington
Seattle, WA
3Division of Environmental Studies
University of California
Davis, CA
4Centre for Ocean Climate Chemistry
Institute of Ocean Sciences
Sidney, BC, Canada
Although the model simulation of pCO2 is satisfactory, this study illustrates the limitations of modeling the chemistry of the upper ocean in one dimension. The slow currents and horizontally homogeneous sea surface in the subarctic Pacific make Papa one of the best available candidates for modeling in 1-D. In spite of this, a 1-D formulation is inadequate to address the chemistry of the halocline (a crucial lower boundary condition to the mixed layer) and does not constrain the transport of the nutrients by wind-driven currents in the mixed layer. We conclude that further progress in modeling the upper ocean nutrient and carbon cycles will require simulation in three dimensions.
A parallel goal for understanding ocean uptake of carbon will be to construct a global scale model of the processes that control sea surface pCO2, that is consistent with our understanding of the physical, chemical, and biological dynamics of the upper ocean. The first step toward a global-scale model is to replicate physical and chemical observations in a few locations where adequate data exist to constrain and test the model formulation.
We describe a physical and chemical model of the upper ocean based on data from Weathership Station Papa in the western subarctic Pacific (50deg.N, 145deg.W). A time-series data set of physical and meteorological observations spanning nearly 30 years exists from this site, recorded by the Canadian weathership program (Tabata (1961)). Because of the regional homogeneity and strong wind forcing of the mixed layer dynamics of this region, many models of upper ocean mixing have been developed or tested based on data from Papa (Denman et al. (1973): Davis et al. (1981): Davis et al. (1981): Martin (1985): Martin (1986): Gaspar (1988): Gaspar et al. (1990): Thomas et al. (1990): Garcon et al. (1992): Thomas et al. (1993)). Chemical studies at the location include a 5-year dataset of pCO2 measurements taken by Wong et al. (1990), associated NO3 measurements (presented in this work), and data taken in conjunction with the SUPER program (Miller et al. (1991); Emerson et al. (1991)).
The model departs from previously published 1-D models of gas chemistry at Papa (Gaspar et al. (1990), Thomas et al. (1990), Garcon et al. (1992), and Thomas et al. (1993)) by including (1) nitrate chemistry, a powerful constraint between mixing and biological production, (2) argon chemistry, which is compared with recent results from the SUPER program (Emerson et al. (1991)), (3) fresh water fluxes and salinity, which constrain the rate of eddy mixing in the halocline. We also include the effects of (4) the Ekman vertical velocity, and (5) we run the model to a long-term steady state (> 30 years), which demonstrates the requirement for isopycnal ventilation of the halocline oxygen and nutrient concentrations. Horizontal advection of nutrients is considered and found to be potentially significant, even conceivably dominant, to the nutrient dynamics at Papa; however, a 1-D formulation is inadequate to address this possibility and we neglect sea surface advection for this study.
An overview of the model formulation is presented in the first section. The calculation of the fluxes of heat and gases at the air-sea boundary is presented in the second section. Since chemical tracers help to constrain the rate of physical mixing, the representation of the chemical processes of new production and respiration is explained in the third section, before an explanation of the physical processes of mixing, upwelling, and horizontal processes in the fourth section. The model behavior is compared with observations in the last section.
Our model is formulated in one dimension, from the sea surface to 198 meters depth, on a vertical grid of 3 meters. The heat fluxes are imposed at the surface based on calculations from the weathership data. The model timestep is 3 hours (based on the timing of meteorological observations). The dynamics of the mixed layer are predicted using the shear instability model of Price et al. (1986). The bottom boundary conditions for T, S, and chemical solutes (O2, Ar, alkalinity, and total CO2) are the observed concentrations at depth of 198 meters; the values above this depth to the sea surface are determined by the model.
Removal of organic carbon and nutrients from the euphotic zone is imposed based on measurements of new production made in conjunction with the SUPER program (Emerson et al. (1991)). For comparison, we also show the model results for an abiotic ocean. Two parameterizations of gas exchange are compared: the relation of Liss et al. (1986), which was found to be consistent with the radon data of Emerson et al. (1991), and the relation of Wanninkhof (1992). We also show results with and without gas transfer by bubble injection as parameterized by Spitzer et al. (1989). A summary of model runs is presented in Table 1.
The model is ground-truthed to the observed T, S, and O2 time series data from the weathership program and dissolved O2 and Ar data from Emerson et al. (1991). The model sea surface pCO2 is compared to the pCO2 data of Wong et al. (1990) for the time period of 1973-1978.
Sensible Heat Flux. The conductive, or sensible, heat flux was calculated using
where [[rho]]air is the density of air, Cp(air) is the specific heat of air at constant pressure, CH is a transfer coefficient equal to 1.1 . 10-3 or 0.8 . 10-3 in unstable or stable atmospheric conditions, respectively (Gill (1982)), U10 is the wind speed (m s-1) at a height of 10 m, SST is the sea surface temperature (deg.C) and Tair is the air temperature. The air temperature and wind speed were taken from the weathership record. SST is predicted by the model.
Latent Heat Flux. The evaporation rate was calculated using the relation
where the transfer coefficient, CE, is taken to be 1.5.10-3 (Gill (1982)), and Qaw and Qair are the specific humidity of the air at the water surface and in the bulk boundary layer, respectively. The humidity at the air/sea interface is taken to be the saturation value at the temperature of the sea surface. The overlying air humidity can in theory be calculated from wet and dry bulb temperature values from the weathership data, but this data implied relative humidity values of 30-50%, which is lower than expected for the marine boundary layer. We used the relative humidity measurements made as part of the pCO2 monitoring effort of Wong et al. (1990). They found relative humidity values of 84 +/- 9 (2[[sigma]])%, with little seasonal cycle. The annual average evaporation rate by this calculation is 0.66 m yr-1, consistent with the estimate of Tabata (1961). The latent heat flux was calculated from the predicted evaporation rate using
QL (W m-2) = (2500.8 + 2.3 SST) Evap
where the term in parentheses on the right hand side is the latent heat of vaporization of water in joule kg-1 (Gill (1982)).
Fresh Water Balance. A calculation related to the latent heat flux is the calculation of the fresh water balance (evaporation - precipitation). Precipitation data were recorded by the weatherships using rain gauges (Knox (1991)). The measured annual precipitation rate is presented in Figure 2. According to the data set, the precipitation rate during the period of interest (1973-1978) was lower than the calculated evaporation rate (solid line, Figure 2). The rainfall record from Victoria, B.C. is generally consistent with the O.W.S. Papa data set for the time period 1954-1980, but the agreement for the modeled period, 1973-1978 is not as good (Figure 2, replotted from Knox (1991). The excess of evaporation over precipitation is impossible to reconcile with preservation of the salinity minimum at the surface. Horizontal transport by wind-driven currents in the mixed layer may represent a further source of salinity to Papa (Tabata (1961); also see below). However, our assessment of the advective salinity flux is based on a long-term climatology, which may not be representative of the period 1973-1979. Royer (1982) concludes that rainfall and runoff at the south coast of Alaska were higher in the middle 1970's than the 1960's by ~20%.
The values of the upwelling rate and the excess freshwater input at the sea surface combine to determine the difference in salinity between the surface and depth. Using an annual average upwelling rate of 27 m yr-1 from Han et al. (1981) (see below), assuming that the upwelling water comes from 200 meters (an assumption that is not at all clear), and accepting that the evaporation rate is known, we require a precipitation rate of 1.7 m yr-1, an increase from the weathership record by a factor of 2.6 during the period from 1973 to 1978. This is consistent with the atlases of Peixoto et al. (1983)and Baumgartner et al. (1975). The fresh water balance presented by Tabata (1961) indicated a slight excess of precipitation over evaporation during the years 1957-1959 (a period of higher rainfall than 1973-1979). Previous models of upper ocean chemistry at Papa (Gaspar et al. (1990), Thomas et al. (1990), Garcon et al. (1992), and Thomas et al. (1993)) avoid the issue by using a time-independent climatological salinity profile, and neglecting the sea surface fresh water fluxes.
Radiation. The incoming radiative heat flux was directly measured as part of the weathership program. These data were found to give a significantly better simulation of the annual cycle of upper ocean heat content (see below) than did estimated solar heat fluxes based on cloud data and the empirical formula of Dobson et al. (1988). The out-going longwave radiative heat flux was calculated as
where [[epsilon]] is the emmisivity (= 0.97), [[sigma]] is the Boltzmann constant (equal to 5.6699 .10-8), Vp air is the vapor pressure of the boundary layer (g m-3), and Cld is the recorded weathership cloud fraction in octas (Gill (1982)).
Uncertainties in Heat Fluxes. Estimation of air-sea heat fluxes involves uncertainties that are typically large relative to the precision desired. The situation is at an optimum in the present study, because of the existence of the weathership data. The uncertainties in the total heat fluxes calculated by bulk formulae have been estimated to be on the order of 25 W m-2 (Weare et al. (1981); Seager et al. (1988); Davis et al. (1981)); roughly 20% of the net fluxes at Papa.
By comparing the heat content of the upper water column between model and observations, we affirm that the calculated heat fluxes are close to the correct values (Figure 3). It must be remembered that most of the calculations depend on the model SST, which is determined in part by mixing, described below.
Gas Exchange. The kinetics of the gas transfer from the atmosphere to the dissolved state in the ocean is typically parameterized as a function of the wind speed. Parameterizations of gas exchange rate have been published by Liss and Merlivat (1986), Tans et al. (1990), and summarized more recently by Wanninkhof (1992). At high wind speeds typical of Papa, the Liss and Merlivat formulation is 30-50% lower than either of the other two formulations. Emerson et al. (1991) measured profiles of dissolved radon and radium at Papa and found that their calculated gas exchange rates were indistinguishable from the Liss and Merlivat formulation, and significantly lower than predicted by either the Tans or Wanninkhof formulations (Figure 4). Accordingly, we use the Liss and Merlivat formulation, which is
where Kw is in units of piston velocity, cm hr-1. Sc is the Schmidt number, calculated as
Sc = [[nu]] D.
where [[nu]] is the viscosity of water and D is a molecular diffusion coefficient. Using this gas exchange relation, we calculated the air-sea flux of oxygen, for example, using
The relative saturations of oxygen, argon, and CO2 were calculated using their solubilities (Weiss (1970); Weiss (1974)) based on the SST, and the atmospheric partial pressure based on the observed barometric pressure from the weathership meteorological data.
A further complication to modeling dissolved gases in the upper ocean is the injection of atmospheric gases into the ocean by collapsing bubbles. This process should lead to a steady-state supersaturation of all of the atmospheric gases. The effect of bubble injection can be estimated by measuring the relative supersaturations of N2 and Ar (which have very different solubilities). At Papa, Emerson et al. (1991), concluded that the observed gas injection rate was consistent with the formulation of Spitzer et al. (1989), which we used here:
where the flux, F, is given in m3 gas at STP m-2 s-1, [[alpha]]i is the mole fraction of the gas in air, [[alpha]]inj is the amplitude of air injection (0.39), ft and fp are the relative fractions of total and partial bubble collapse (0.07 and 0.93), [[Gamma]] is the ratio of the gas diffusivity to that of helium, and [[gamma]] is an exponent of 2.2.
Unfortunately, these measurements were all made in the summer, and the winter new production rates are largely unconstrained by recent data. Measurements of net production by 14C incubations reveal lower wintertime values by a factor of 4 (Parsons et al. (1988)). These data were taken before the development of clean techniques, and were therefore a factor of 2-4 lower than values measured today (Booth et al. (1988)). However, since the early data is the only indicator for winter values, we assume that winter organic carbon export rates are a factor of 4 lower than summer values.
The flux of sinking particles measured in sediment traps is generally found to decrease as a function of trap depth (Martin et al. (1987)), as:
where Foc(z) is the sinking flux of organic carbon at depth z. This empirical observation is consistent with Welschmeyer (unpublished) sediment trap data from Papa. Although some of the sinking flux measured near the sea surface may be an artifact of biological activity ("swimmers": Lee et al. (1988)) we assume that the decrease with depth represents in situ degradation of sinking organic matter. Since we take the depth of the euphotic zone to be 50 meters, we also assume that the sinking flux of carbon is attenuated in the region from 50-100 m using the same relation. The implication of this relation is that 50% of the carbon flux exported from the euphotic zone is remineralized before it sinks to 100 m (and hence the redissolved nutrients would be available for re-entrainment in the winter), and only 25% of the total carbon export sinks past 200 m depth. Respiration should represent a source for CO2 and NO3, and a sink for O2. However, we find that the respiration sink for O2 in the subsurface waters (100-200 m) must be balanced by isopycnal ventilation, which is neglected in a 1-D model formulation (see below).
The carbonate buffer system in the model was treated using the full pH equilibrium system of CO3, HCO3, and CO2. Profiles of alkalinity and total CO2 were maintained at all depths, along with the other quantities (T, S, O2, Ar, and NO3). For a given time step, the pCO2 in the surface box was calculated, from which we calculated the air-sea flux of CO2. Using a bottom boundary condition of a constant observed concentration (alkalinity and total CO2 from Cline et al. (1985), the system was run to steady state.
where C is the scalar quantity, h is the depth of the mixed layer, wh is the vertical velocity at the base of the mixed layer, Cw is the scalar quantity in the mixed layer if wh is downward and the value immediately below the mixed layer if wh is upward [[kappa]] is an eddy diffusion coefficient, ([[partialdiff]]C/[[partialdiff]]z)h is the scalar quantity depth gradient at depth h (just below the mixed layer), S is the surface flux (for examples atmospheric heat exchange, gas exchange; disscussed below), and B is a source or sink internal to the mixed layer (the only example of this is internal heating by absorbtion of solar radiation).
Three general families of mixed layer models have been developed to predict the time evolution of the mixed layer depth given atmospheric forcing (wind velocity and air-sea heat fluxes). Although the assumed physics controlling mixing differ among the three types, all are reasonably successful at predicting mixed layer depth fluctuations. The three families are the "integrated turbulent kinetic energy" models (Kraus et al. (1967), the "turbulence closure" models (e.g. Mellor et al. (1974)), and the "shear instability" models (e.g. Price et al. (1986)). The relative merits and performances of the three models have been compared elsewhere (Martin (1985); Martin (1986); Archer (1990); Thomas et al. (1993)). We have chosen the Price et al. (1986) shear instability model because it is simple, as accurate as any of them, and because the empirical parameters required to predict mixing (see below) seem to be constant over a range of environmental conditions (Price et al. (1978)).
The model is based on the premise that energy for mixing is derived chiefly from shear at the mixed layer base. Wind energy transfers momentum into the upper ocean in the form of mean current, which rotates inertially. Mixing is predicted using a non-dimensional parameter called the bulk Richardson number
where g is the acceleration due to gravity, [[Delta]][[rho]] and [[Delta]]u are the difference in density and velocity between the mixed layer and the layer immediately below it (the vertical resolution in the model was 3 m), and h is the mixed layer thickness. When Rib is less than a "critical" value of 0.65, the properties (T, solutes, and momentum) of the mixed layer are averaged downward (completely mixed) until the new Rib becomes higher than the critical value.
Two other mechanisms for mixing are included in the original Price et al. (1986) model. Convective mixing relieves the condition of static instability that results from negative buoyancy flux at the sea surface; this aspect of the model is retained. Price et al. (1986) also included a shear dependent mechanism for partial mixing beneath the homogeneous mixed layer; this part of the formulation was replaced, as described below.
A further uncertainty is the origin of the upwelling water; upwelling from 200 meters, as is assumed here, will carry different chemical load than if the water is delivered isopycnally from a shallower depth in the halocline.
For salinity, there appears to be an along flow gradient of perhaps 0.25 salinity units on a distance scale of 1000 km (to the bend in the flow at 45N, 155W). If we assume that the annual average mixed layer is 40 meters deep, we can calculate that the advective flux of salinity is equivalent to 1.5 meters of evaporation from the sea surface. Even without this term, there appears to be greater flux of salt into the model than can be balanced by precipitation, during these years of low rainfall (see Fresh Water Balance section, above). Including an advective source of salt will only make the problem worse. For dissolved phosphate, we calculate that a concentration differential of 0.15 uM PO4 from 45N, 155W would be sufficient to supply all of the nutrient required to balance the annual new production, with no contribution from upwelling nutrient-rich halocline water at all.
We feel that the potential importance of horizontal advective fluxes of salt and nutrients is a manifestation of the limitations of one-dimensional models for simulating upper ocean chemistry. Station Papa is an optimum candidate for a one-dimensional formulation, since the currents are slow and the chemistry horizontally homogeneous relative to most other parts of the world's oceans. Even here the limitations of a 1-D model are apparent. Since the model is unable to constrain the advective sources and sinks in the surface layer, we feel that to attempt to include these affects would sacrifice the integrity of the model to dominance by tunable parameters. For example, modeling the dissolved gases would require including or parameterizing the time history of sea surface temperature. Accordingly, we neglect these fluxes, and continue with our one-dimensional formulation, leaving a full 3-D simulation for future work. Previous models of the chemistry of the upper ocean have simulated the dynamics of oxygen and carbon uncoupled to the supply of nutrients (the models did not include scaler nitrate or phosphate fields), and so the advective nutrient fluxes were neglected without comment Thomas et al. (1990), Gaspar et al. (1990), Thomas et al. (1993).
Vertical Mixing. To the extent that mixing in the waters below the mixed layer is dominated by vertical processes, the equation governing the distribution of a scaler quantity (heat, salinity, and dissolved tracers) is:
The rate of eddy mixing ([[kappa]]h) in the halocline is not tightly constrained by observational data. The original Price et al. (1986) formulation included a mechanism for partial mixing in the stratified region based on the gradient Richardson number. In practice, this mechanism in the model only mixes a few meters below the bottom of the completely mixed layer, and has little effect on turbulent transport of solutes through the halocline and entrainment into the mixed layer Large et al. (1986). Mixing below the shear zone at the mixed layer base is thought to be driven by internal wave breaking, and is not addressed by the Price et al. (1986) model. Denman et al. (1988) estimated that Kv, the eddy diffusivity, was in the range of 0.03 - 0.3 cm2 s-1 in the seasonal thermocline at station Papa based on microstructure velocity measurements. Miller et al. (1991) estimated halocline diffusivity of 0.6 cm2 sec-1. Based on the evolution of profiles of temperature with depth taken from free-drifting thermister chains, Large et al. (1986) find that periods of intense mixing through a depth of ~1.5 mixed layer depths often accompany high-wind events. Their data were best fit using a value of 4 cm2 s-1 during storms, and a "background" value of 0.2 cm s-1 during periods of calm.
Our choice of eddy mixing coefficients is motivated by two constraints: the shape of the salinity profile with depth, and the depth of the winter mixed layer. We find that, although the heat fluxes appear to be correct (Figure 3), the winter mixed layer depth is underestimated by the model. The discrepancy is helped by incorporating storm-driven events of high mixing, such as observed by Large et al. (1986). We define storms as periods of winds exceeding 28 knots for a duration 18 hours or longer. There are a total of 139 such events during the 6 year record between 1973 and 1978. During storms, we assign Kv a value of 4 cm2 s-1, following Large et al. (1986).
However, if we assume that the effect of isopycnal exchange on the salinity structure is second-order, then the shape of the salinity profile becomes another constraint on Kv (Figure 6). In order to reproduce the observed curvature of the salinity profile, we require a "background" Kv of ~0.2 cm2 s-1, during periods that are not storms. This constraint would have been missed, had we neglected to predict salinity and fresh water fluxes at the sea surface.
Isopycnal Ventilation. Although the salinity structure of the halocline suggests that vertical processes are important to the distribution of salt, we find that this formulation alone is inadequate to explain the distribution of oxygen and nutrients. Recall our assumption that the redissolution of sinking organic particles represent a sink for O2 in the halocline, whose magnitude can be roughly estimated as the divergence of sediment trap fluxes with depth. We find that vertical diffusion by itself is inadequate to replenish the oxygen consumption in the halocline. After ~ 30 years of model simulation, the halocline becomes anoxic (Figure 7). Given an observed steady-state vertical gradient of oxygen of 150 uM / 100 m, and a consumption rate subsurface roughly equal to the generation rate in the euphotic zone (here taken to be of order 10 mM m-2 d-2), a diffusion coefficient of approximately 10-4 m2 sec-1 is required to balance oxygen consumption by vertical eddy mixing. The eddy diffusion coeffient appropriate here was constrained to be lower than that by the shape of the salinity profiles (see section above).
We assume that the oxygen source is ventilation along isopycnal surfaces. The density of the water at 200 meters is approximately 26.8, the densest water that directly outcrops in winter (Reid (1965)). Isopycnal ventilation is consistent with freon data (Warner (1988)) and the time evolution of the tritium distribution (VanScoy et al. (1991)). Jenkins (1982) conceptually balanced in situ oxygen consumption with isopycnal ventilation, rather than with local processes, to calculate oxygen consumption rates. However, previous models of the nutrient and oxygen cycles in the upper ocean have neglected isopycnal processes (Peng et al. (1987), Musgrave et al. (1988), Thomas et al. (1990), Garcon et al. (1992), and Thomas et al. (1993)). These efforts apparently evaded the requirement for isopycnal subsurface ventilation by using a higher eddy diffusion coefficient (see section above), a lower rate of production, or by spinning up the simulations for only one or a few years. It must be noted, however, that the winter mixed layer in the North Pacific is much shallower than in the Atlantic, due to the salinity stratification; isopycnal ventilation may play less of a role in the upper ocean nutrient dynamics therefore in the Atlantic.
The established.gifure of ventilation of the North Pacific was given by Reid (1965). According to his analysis, waters of this depth are ventilated from the Northeast Pacific near the Sea of Okhotsk. The required rate of ventilation can be estimated to first order by assuming a balance between removal and supply of O2 in the halocline waters. Based on the empirical organic carbon redissolution relation of Martin et al. (1987) (see above), we estimate that something like 25% of the flux past 50 meters is redissolved between the depths of 100 and 200 meters. Using an annual average of ~90 mg C m-2 d-1 carbon export from the euphotic zone, we calculate that the average oxygen demand for the waters between 100-200 m depth is 0.009 mol O2 m-3 y-1. The diffusional and advectional terms in the oxygen balance are small in this region, where the observed gradient is small. Therefore, this O2 consumption term must be largely balanced by isopycnal ventilation. If we formulate ventilation as the simple exchange of water between oxygen-saturated source water and halocline Papa water (with an oxygen concentration contrast of 100 uM), then we can calculate that the halocline waters at Papa must be ventilated on a time scale of ~10 years. This is in the same range as the estimate of VanScoy et al. (1991), based on the time evolution of the tritium distribution, of "<14 years".
We have demonstrated the inadequacy of modeling the dissolved oxygen chemistry of the thermocline in a one-dimensional sense, even in the 100 meters immediately below the wintertime mixed layer. Unfortunately, the process of isopycnal ventilation is not easy to constrain in a one-dimensional formulation. We treat dissolved oxygen (and, analogously, nitrate and carbon) as if the local sink by respiration is exactly balanced by isopycnal ventilation. In other words, the respiration oxygen sink is simply neglected, as if sinking organic particles were refractory to respiration. Carbon and nutrients removed from the euphotic zone by production are not regenerated at depth, but simply removed from the simulation. Unfortunately, since the processes of respiration and ventilation are not addressed in detail, the model is unable to constrain independently the flux of nitrate into the euphotic zone. The halocline part of the model has become a support or boundary condition for the mixed layer and seasonal thermocline, rather than predictive by itself, which will have to await a 3-dimensional formulation. Although we have little faith in the mechanics of the model in the halocline, we find that the supply rate of nitrate to the euphotic zone balances the observed new production rate, using this model configuration (note: a further problem is the advective source of nitrate to Papa by wind-driven currents; see above).
However, since the temperature and salinity structure resemble the expected results of a vertically forced system (Tully et al. (1960)), we assume that isopycnal advection plays only a secondary part in determining T and S. The diffusion coefficient and upwelling rate required to maintain salinity profiles that resemble reality are not outside the bounds of the expected values, as outlined above.
The simplest comparison is to look at the oxygen concentration by itself. An analysis of this type was used by Musgrave et al. (1988), and by Thomas et al. (1990) and Emerson (1987), to estimate the new production rate at Bermuda and at Papa, respectively. We compare the results of the Standard, RW, No Bubbles, and Dead model simulations (see Table 1) with the observations from the weather ship data in Figure 11 (a). It appears that the sea surface oxygen concentrations of the model are insensitive to gas exchange rates, bubble injection, and biological activity relative to the variations in the measured data.
The effects of temperature and biology can be separated by comparing concentrations of oxygen and argon simultaneously. These gases have very similar temperature dependence of their relative saturation values;
where [[Delta]]O2 is the concentration of oxygen relative to the saturation value, etc. (Benson (1964)). A similar relation can be found for pressure (Benson (1964))
F([[partialdiff]][[Delta]]O2, [[partialdiff]]P) = 1.04 F([[partialdiff]][[Delta]]Ar, [[partialdiff]]P).
where changes in pressure can be caused by changes in atmospheric pressure or by injection of bubbles into the water column. Thus, to a large extent, argon acts as an abiotic analog for oxygen, so that the effects of physical forcing can be removed from the dissolved oxygen signal by subtraction of the argon signal. A plot of [[Delta]]O2 - [[Delta]]Ar is presented in Figure 11 (b). For the "Dead" model run, the value of [[Delta]]O2 - [[Delta]]Ar is close to zero during summer, as expected if the physical effects on the two gas concentrations are indeed similar. In the winter, mixing from below of oxygen deficit (the bottom boundary condition is the same for all three model runs) lowers the oxygen concentration relative to argon (which has a bottom boundary condition of the saturation value). Data from Emerson et al. (1991) from the SUPER program is also shown on Figure 11 (b). Although the data reflect large natural variability, Emerson's conclusion is that the new production rate at Papa is on the order of 150 mg C m-2 day-1: the same as used here for the Standard, RW, and No Bubbles runs (Table 3). The oxygen concentration in the water column from the Standard model are compared with weathership data in Figure 11 (c).
pCO2. The model sea surface pCO2 time series is compared with data of Wong et al. (1990) in Figure 13. The atmospheric value was generated by taking a composite year from the atmospheric values in the Wong dataset. Comparison of the model results (solid line) with data (+) reveals that (1) the model appears to get the annual mean value correct (see Table 4), relative to the range of observed pCO2 in the world ocean, which is ~150 -550 ppm (Tans et al. (1990)). (2) The annual cycle appears to be simulated by the model, with values of ~350 ppm in late summer decreasing to ~300 ppm in winter. (3) Short-term fluctuations in the data do not appear to be captured by the model pCO2 curve, which appears to be steadier than the data.
We believe that the fast fluctuations in the observed pCO2 are tied to fluctuations in the observed concentrations of nutrients in the mixed layer. To verify this possibility, we ran the model in a "NO3 Restoring" mode, where the new production value was calculated such that the model nitrate sea surface nitrate concentration corresponds to observations (Figure 10). During periods of missing data, the model new production was taken to be the specified values of 150 mg C m-2 day-1 in summer and 37 mg C m-2 day-1 in winter, as in the Standard simulation. Although the variation in the observed NO3 is more likely to be caused by patchiness in the NO3 distribution than by rapid changes in uptake and removal, the effects on the dissolved carbon concentration should be similar. The resulting model surface NO3 time series is presented as the dashed line in Figure 10. The model pCO2 (dotted line in Figure 13) exhibits the same fast fluctuations as observed in the data, confirming that fluctuations (patchiness) of the nutrient concentration could be the cause of the spikey character of the pCO2 data relative to the standard model.
Two alternate models of biological activity are contrasted with the Standard model in Figure 14 (a). The top curve is was generated by the Dead model (see Table 1 for a description of model experiments). The value of the sea surface pCO2 from this model is higher than in the standard model by approximately 60 ppm. The lower curve is the result of the Bloom model, in which the biological system removes the sea surface nutrients to zero values, with a time constant of two weeks, when the euphotic zone is shallower than 50 meters. The average pCO2 of the bloom model is similar to the Standard model, but with negative pCO2 spikes of 20-40 uatm in summer. The air - sea fluxes of CO2 for the three models are presented in Figure 14 (b) and Table 4. While the Standard model, designed to replicate station Papa observations, has near zero net flux of carbon across the air-sea interface, the bloom model takes up carbon, and the dead model is a source of carbon to the atmosphere. The simulations demonstrate the importance of the upper ocean biological system on determining the pCO2 of the surface water.
Dissolved oxygen is consumed below the mixed layer by respiration of sinking particles (as inferred from the decrease in Welschmeyer's (unpublished) sediment trap data with depth). Maintenance of an oxic halocline requires isopycnal ventilation, with a required time constant of ~11 years, consistent with previous estimates based on tritium and freon data. Gas exchange at the sea surface was parameterized using the wind-speed dependence of Liss et al. (1986). The gas exchange formulation was ground-truthed against dissolved oxygen and argon data. The time series of dissolved oxygen was slightly closer to saturation than the weathership hydrographic data. Oxygen/argon relative saturation values were consistent with measurements of Emerson et al. (1991). Wind-driven currents at the sea surface represent a potential source of dissolved nutrients that is equal in magnitude to the flux by upwelling. We conclude that this possibility, and the problem of isopycnal ventilation, demonstrate the limits of the utility of a one-dimensional model for the chemistry of the upper ocean.
The model sea surface pCO2 appeared to capture the annual mean and seasonal cycle of the observed data of Wong et al. (1990). The model pCO2 was much smoother than the data on time frames of a few weeks, however; by running the model in a "nutrient restoring" mode (where model surface NO3 concentrations are tied to observations) we are able to replicate the observed pCO2 fluctuations. We assume that the observed fluctuations in NO3 and pCO2 reflect patchiness in the upper ocean. Changes in the upper ocean biological model, to a "bloom" system typical of the North Atlantic, or a "dead" ocean, have a significant impact on the steady state pCO2 of the model sea surface.
ACKNOWLEDGMENTS. This paper benefited from discussions with Jim Price, Taro Takahashi, Bruce Frost, Gilles Reverdin, and James Bishop, and was supported by the Lamont Doherty Geological Observatory, the Department of Energy contract no. DE-AC06-76RLO1830, the G. Unger Vettersen Foundation, and selected member companies of the International Petroleum Industry Environmental Conservation Association.
Archer, D. E., (1990): Modeling pCO2 in the upper ocean: A review of relevant physical, chemical, and biological processes. Dept. of Energy Tech. Rept. TRO50.
Baumgartner, A. and E. Reichel (1975): The world water balance: Mean annual global, continental and maritime precipitation, evaporation and run-off. Elsevier, Amsterdam.
Benson, B. B. (1964): Some thoughts of gases dissolved in the oceans, In D. R. Schink and J. T. Corliss (Eds.), Proc. Symp. on Mar. Geochem., Univ. Rhode Island,
Booth, B. C., J. Lewin and C. J. Lorenzen (1988): Spring and summer growth rates of subarctic Pacific phytoplankton assemblages determined from carbon uptake and cell volumes estimated using epifluorescence microscopy, Mar. Biol., 98: 287-298.
Broecker, W. S. and T. H. Peng (1982): Tracers in the sea. Eldigio Press, Palisades, NY.
Cline, J. D., R.A.Feely, K. Kelly-Hansen, J. F. Gendron, D. P. Wisegarver and C. T. Chen, (1985): Current inventory of anthropogenic carbon dioxide in the North Pacific subarctic Gyre. NOAA Technical Memorandum ERL PMEL-60.
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 Res., 28A(12): 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 Res., 28A(12): 1453-1475.
Denman, K. L. and A. E. Gargett (1988): Multiple thermoclines are barriers to vertical exchange in the subarctic Pacific during SUPER, May 1984, J. Mar. Res., 46: 77-103.
Denman, K. L. and M. Miyake (1973): Upper ocean modification at Ocean Station Papa: observations and simulation, J. Phys. Ocean., 3: 185-196.
Dobson, F. W. and S. D. Smith (1988): Bulk models of solar ratiation at sea, Q. J. R. Meteorol. Soc., 114: 165-182.
Emerson, S. (1987): Seasonal oxygen cycles and biological new production in surface waters of the subarctic Pacific Ocean, J. Geophys. Res., 92(C6): 6535-6544.
Emerson, S., P. Quay, C. Stump, D. Wilbur and M. Knox (1991): O2, Ar, N2, and 222Rn in surface waters of the subarctic ocean: Net biological O2 production, Global Biogeochem. Cycles, 5: 49-71.
Fabry, V. J. (1989): Aragonite production by pteropod molluscs in the subarctic Pacific, Deep Sea Res., 36: 1735-1751.
Feely, R. A., R. H. Byrne, J. G. Acker, P. R. Betzer, C. T. A. Chen, J. F. Gendron and M. F. Lamb (1988): Winter -- summer variations of calcite and aragonite saturation in the Northeast Pacific, Mar. Chem., 25: 227-241.
Garcon, V. C., F. Thomas, C. S. Wong and J. F. Minster (1992): Gaining insight into the seasonal variablility of CO2 at ocean station P using an upper ocean model, Deep Sea Res., 39: 921-938.
Gaspar, P. (1988): Modeling the seasonal cycle of the upper ocean, J. Phys. Ocean., 18: 161-180.
Gaspar, P., J.-C. Andre and J.-M. Lefevre (1990): The determination of the latent and sensible heat fluxes at the sea surface viewed as an inverse problem, J. Geophys. Res., 95: 16169-178.
Gaspar, P., Y. Gregoris and J. M. Lefevre (1990): A simple eddy kinetic energy model for simulations of the oceanic vertical mixing: tests at station papa and the long-term upper ocean study site, J. Geophys. Res., 95: 16179-193.
Gill, A. E. (1982): Atmosphere-Ocean Dynamics. Acedemic Press, Orlando.
Han, Y.-J. and S.-W. Lee (1981): A new analysis of monthly mean wind stress over the global ocean. Oregon State University, Corvallis, Oregon.
Jenkins, W. J. (1982): Oxygen utilization rates in North Atlantic subtropical gyre and primary production in oligotrophic systems, Nature, 300: 246-248.
Knox, J. L. (1991): An assessment of the 27-year record of measured precipitation at ocean weather station "P" in the northeast Pacific Ocean, Clim. Bull. Canadian Met. Ocean. Soc, submitted:
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.
Large, W. G., J. C. McWilliams and P. P. Niiler (1986): Upper ocean thermal response to strong autumnal forcing the northeast Pacific, J. Phys. Ocean, 16: 1524-1550.
Lee, C., S. G. Wakeham and J. I. Hedges (1988). The measurement of ocean particle flux -- are "swimmers" a problem? Oceanography,
Liss, P. S. and L. Merlivat (1986): Air-sea gas exchange rates: Introduction and synthesis, In Baut-Menard (Eds.), The role of air-sea exchange in geochemical cycling, D. Reidel,
Martin, J. H., G. A. Knauer, D. M. Karl and W. M. Broenkow (1987): VERTEX: carbon cycling in the northeast Pacific, Deep Sea Res., 34(2): 267-285.
Martin, P. J. (1985): Simulation of the mixed layer at OWS November and Papa with several models, J. Geophys. Res., 90(C1): 903-916.
Martin, P. J. (1986): Testing and comparison of several mixed-layer models, Naval Ocean Research and Development Activity Tech Rep, 143:
McNally (1981): Satellite-tracked drift buoy observations of the near-surface flux in the eastern mid-latitude north Pacific, J. Geophys. Res., 86: 8022-8030.
McNally, G. J., W. C. Patzert, A. D. J. Kirwan and A. C. Vastano (1983): The near-surface circulation of the north Pacific using satellite tracked drifting buoys, J. Geophys. Res., 88: 7505-7518.
Mellor, G. L. and T. Yamada (1974): A hierarchy of turbulence closure models for planetary boundary layers, J. Atmos. Sci., 31: 1791-1806.
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, Limnol. Oceanogr., 36: 1600-1615.
Musgrave, D. L., J. Chou and W. J. Jenkins (1988): Application of a model of upper-ocean physics for studying seasonal cycles of oxygen, J. Geophys. Res., 93: 15679-15700.
Parsons, T. R. and C. M. Lalli (1988): Comparative oceanic ecology of the plankton communities of the subarctic Atlantic and Pacific Oceans, Oceanogr. Mar. Biol. Annu. Rev., 26: 317-359.
Peixoto, J. P. and A. H. Oort (1983): The atmospheric branch of the hydrological cycle and climate, In A. S.-P. e. al. (Eds.), Variations of the global water budget, Reidel, London.
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.
Price, J. F., C. N. K. Mooers and L. Van (1978): Observations and simulation of storm-induced mixed-layer deepening, J. Phys. Ocean., 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, J. Geophys. Res., 91(C7): 8411-8427.
Reid, J. (1965): Intermediate waters of the Pacific Ocean. Johns Hopkins Press, Baltimore.
Royer, T. C. (1982): Coastal fresh water discharge in the Northeast Pacific, J. Geophys. Res., 87: 2017-2021.
Seager, R., S. E. Zebiak and M. A. Cane (1988): A model of the tropical pacific sea surface temperature climatoloty, J. Geophys. Res., 93 (C2): 1265-1280.
Spitzer, W. S. and W. J. Jenkins (1989): Rates of vertical mixing, gas exchange, and new production: estimates from seasonal gas cycles in the upper ocean near Bermuda, J. Mar. Res., 47: 169-196.
Tabata, S. (1961): Temporal changes of salinity, temperature, and dissolved oxygen content of the water at station "P" in the Northeast Pacific Ocean, and some of their determining factors, J. Fish. Res. Bd. Canada, 18(6): 1073-1124.
Tans, P. P., I. Y. Fung and T. Takahashi (1990): Observational constraints on the global atmospheric carbon dioxide budget, Science, 247: 1431-1438.
Thomas, F., V. Garcon and J.-F. Minster (1990): Modeling the seasonal cycle of dissolved O2 in the upper ocean at Ocean Weather Station P., Deep Sea Res., 37: 463-491.
Thomas, F., J. F. Minster, P. Gaspar and Y. Gregoris (1993): Comparing the behaviour of two ocean surface models in simulating dissolved O2 concentration at O.W.S. P., Deep Sea Res., 40: 395-408.
Tully, J. P. and F. G. Barber (1960): An estuarine analogy in the sub-arctic Pacific Ocean, J. Fish. Res. Board Canada, 17: 91-112.
VanScoy, K. A., R. A. Fine and H. G. Ostlund (1991): Two decades of mixing tritium in to the north Pacific Ocean, Deep Sea Res., 38: S191-S220.
Wanninkhof, R. (1992): Relationship between wind speed and gas exchange over the ocean, J. Geophys. Res., 97: 7373-7382.
Warner, J. M. (1988): Chlorofluoromethanes F-11 and F-12: their solubilities in water and seawater and studies of the distributions in the South Atlantic and North Pacific Oceans. Univ. Cal. San Diego,
Weare, B. C., P. T. Strub and M. D. Samuel (1981): Annual mean surface heat fluxes in the tropical Pacific Ocean, J. Phys. Ocean., 11: 705-717.
Weiss, R. (1974): Carbon dioxide in water and seawater. The solubility of a non-ideal gas, Mar. Chem., 2: 203-215.
Weiss, R. F. (1970): The solubility of nitrogen, oxygen, and argon in water and seawater, Deep Sea Res., 17: 721-735.
Wong, C. S. and Y.-H. Chan (1990): Temporal variations in the parital pressure and flux of CO2 at ocean station P in the subarctic northeast Pacific Ocean, Tellus, 43B: 206-223.
Table 1. Summary of model runs. LM refers to the gas exchange parameterization of Liss et al. (1986); RW refers to the parameterization of Wanninkhof (1992). SJ refers to the bubble injection parameterization of Spitzer et al. (1989). Values of organic carbon export production are given in units of mg C m-2 day-1, summer/winter, where the summer value is imposed when the mixed layer is shallower than 50 meters depth. The Bloom experiment depleted the available NO3- to a value of 0 with a time constant of 24 days in summer. The NO3 restoring model restored the surface NO3 values to observations with a time constant of 1 day.
Run Gas Exchange Bubble Export Carbon Production
Designation Injection
Standard LM SJ 150/37
RW RW SJ 150/37
No Bubbles LM none 150/37
Dead LM SJ 0/0
Bloom LM SJ [[tau]]=24 day/37
NO3 restoring LM SJ NO3 restored to
observations
Table 2. Measurements of new production made in conjunction with the SUPER program. From Emerson et al. (1991).
Method Flux mg C m-2 d-1 Reference
sediment trap 100-130 Welschmeyer
fluxes1
oxygen gas exchange 130-160 Emerson
NO3 budget2 160-200 Wheeler,
unpublished
1
With additional mesozooplankton migration and particulate buildup terms
added.2[[partialdiff]]NO3/[[partialdiff]]t, with diffusive nitrate flux across the mixed layer base added.
Table 3. Rate of formation of calcium carbonate in the upper ocean at station Papa. From Fabry (1989).
production, mg m-2 d-1 % of total
pteropods 4.4 4-13
foraminifera 8.8 18-30
coccolithophorids 36-72 59-77
total 50-85
Table 4. CO2 fluxes between air and sea, for various configurations of the model. "Data" average [[Delta]]pCO2 was calculated from the measurements of Wong and Chan (1991), binned by julian day and averaged. See Table 1 for description of the "various experiments.
Experiment Average Net CO2 Flux,
[[Delta]]pCO2, mmol m-2 d-1
ppm
Standard -10 -0.95
RW -3 -0.73
Dead 42 1.96
Bloom -18 -1.32
No Bubbles -6 -0.71
Restore 2 -0.11
Data 13