Vertical mixing effects on the phytoplankton bloom - University of ...

JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 111, C03002, doi:10.1029/2005JC002994, 2006

Vertical mixing effects on the phytoplankton bloom in the southeastern Bering Sea midshelf Meibing Jin,1 Clara J. Deal,1 Jia Wang,1 Nori Tanaka,1 and Moto Ikeda2 Received 7 April 2005; revised 14 July 2005; accepted 17 November 2005; published 2 March 2006.

[1] A vertically one-dimensional ecosystem model was developed and applied to the

southeastern Bering Sea middle shelf. The physical model includes a 2.5-level turbulence model. Tidal forcing was introduced to improve representation of tidal mixing in addition to wind stirring and thermal stratification. The simulated currents, thermocline and mixed-layer depth (MLD) agree well with observations. The biological model was adapted from Eslinger et al. (2001) with nine compartments. The onset, magnitude and duration of the spring phytoplankton blooms were realistically reproduced at 12 m, 24 m, and 44 m in the standard run. The spring phytoplankton bloom was dominated by diatoms, and the post blooms by flagellates, which agree with previous studies in the region. The peak phytoplankton biomass reached 8 mmol N m3, or approximately 16 mg Chl m3, comparable to that observed in the PROBES program in 1980 and 1981 (Eslinger and Iverson, 2001). The simulated primary productions were within the range of 60 to 200 g C m2/yr estimated in previous studies. Sensitivity studies revealed distinct effects of tidal mixing, wind stirring, thermal stratification and their impacts on the timing and magnitude of the phytoplankton bloom and the gross and net primary production. Links of MLD with primary production and species were found. If a constant MLD is used in the model, there exists a maximum gross primary production (GPP) at MLD = 24 m. Model results reveals that the predominant phytoplankton species changes from flagellates when MLD < 15 m to diatoms when MLD > 15 m. Citation: Jin, M., C. J. Deal, J. Wang, N. Tanaka, and M. Ikeda (2006), Vertical mixing effects on the phytoplankton bloom in the southeastern Bering Sea midshelf, J. Geophys. Res., 111, C03002, doi:10.1029/2005JC002994.

1. Introduction [2] The southeastern Bering Sea is the most productive marine ecosystem in the United States and one of the most productive in the world [Springer et al., 1996]. Any change to this rich ecosystem that causes a reduction in productivity, change in species composition, or change in the portion of the food web that is usable by mankind will have a severe societal impact [Napp and Hunt, 2001]. Because phytoplankton can fix inorganic carbon, phytoplankton productivity also plays an important role in modifying the carbon dioxide content of the atmosphere and hence the scenarios for global climate change. The motivation of this study is to understand how physical processes influence phytoplankton productivity. [3] The NOAA/PMEL biophysical mooring site M2 (Figure 1), the region of interest for this model application, is located in the southeastern Bering Sea middle shelf domain. The three domains (coastal, middle and outer) indicated in Figure 1 are differentiated by hydrographic 1 International Arctic Research Center, University of Alaska, Fairbanks, Alaska, USA. 2 Graduate School of Environmental Earth Science, Hokkaido University, Sapporo, Japan.

Copyright 2006 by the American Geophysical Union. 0148-0227/06/2005JC002994$09.00

features associated with characteristic bathymetric ranges [Hunt and Stabeno, 2002]. The region’s weak net circulation and small heat advection and diffusion provide an ideal site for vertically 1-D physical-biological model studies. The continuous biophysical mooring observations from 1995 to the present illustrate that the Southeastern Bering Sea shelf is a highly variable system at annual and interannual scales [Overland et al., 2001; Schumacher and Alexander, 1999; Francis et al., 1998]. The primary production responds to the varying physical environments through the interplay of light, advection, stratification, nutrient supply, and vertical mixing [Springer et al., 1996]. Secondary production is in turn strongly dependent on the magnitude, quality and timing of this primary production. Spring phytoplankton blooms are generated when organisms are confined to the euphotic layer by vertical stratification developed in the water by melting ice (early ice-related bloom) or the sun’s warming (openwater spring bloom). Hunt and Stabeno [2002] reviewed the time series of mooring data from 1995 to 2001 and found that retreat of the winter sea ice before mid-March (or failure of ice to be advected into the region) results in an open water bloom in May or June in relatively warm water (3C), and conversely, when ice retreat is delayed until mid-March or later, an ice-associated bloom occurs in cold water (0C) in early spring. The sea ice that reaches the mooring site is annual ice that forms in the lee of the islands

C03002

1 of 20

C03002

JIN ET AL.: MIXING EFFECTS ON PHYTOPLANKTON BLOOM

Figure 1. Topography of southeastern Bering Sea. The NOAA-PMEL biophysical mooring site M2 is marked with an asterisk in the middle domain. and coasts of the Bering Sea and is advected southward by the prevailing winds into the mooring site area in winter. Spring blooms at the Bering Sea ice edge have been observed and reported by Niebauer et al. [1990, 1995] and Stabeno et al. [1995]. Since the ice advection process was not included in the vertically 1-D model, in this study we focused on the physical and biological process during the open-water bloom of the year 2000. [4] Meteorological forcing and vertical mixing processes were found to be critical to the spring plankton bloom onset and progression observed during the Processes and Resources of the Bering Sea shelf (PROBES) program in 1980 and 1981 [Eslinger and Iverson, 2001]. Vertical mixing is typically a combination of wind stirring, tidal mixing and thermal mixing/stratification. The middle domain (depth between 50 m and 100 m) is strongly influenced by tidal mixing, as well as wind- and buoyancy-driven subtidal flows and their interactions. Tidal and wind mixing is insufficient to stir the entire water column in the presence of a positive buoyancy flux from spring to early autumn [Stabeno et al., 2001]. Earlier ecosystem models in the region, such as those developed by Eslinger and Iverson [2001] and Merico et al. [2004], are all 1-D coupled physical-biological models, in which the tidal mixing process was parameterized rather than computed explicitly in those models. Merico et al. [2004] used a two-layer model in which tidal mixing was a constant enhancement to the upper mixing layer. Eslinger and Iverson [2001] and Eslinger et al. [2001] computed wind and thermal mixing effects through a Froude number-based model [Pollard et al., 1973; Thompson, 1976], and tidal mixing was parameterized as a constant 40-m benthic boundary layer and a 0.1 h1 mixing rate in both bottom and surface mixing layers. Tidal mixing mainly contributes to formation of the bottom boundary layer, and it affects biology by bringing up rich nutrients from the bottom through vertical mixing. The magnitude of tidal mixing displays variations from tidal period (semidiurnal to diurnal) to seasonal scales due to the interaction with seasonal stratifications. Thus tidal mixing is

C03002

important from semidiurnal and seasonal timescales. Although tidal mixing can be readily simulated with a 3-D ocean model in the Bering Sea, such as that of Hermann et al. [2002], the tidal current has not been implemented in a vertically 1-D physical-biological model owing to the difficulty of applying tidal forcing in the 1-D model setting. Therefore the influences of tidal mixing on the primary productions have not been reported, although vertically 1-D physical-biological models have also been widely used in various regions for ecosystem studies [e.g., Kuhn and Radach, 1997; Chai et al., 2002; Fujii et al., 2002]. [5] The continuous, vertically fine-scale biophysical mooring data at M2 make it possible to develop and validate a 1-D coupled physical-biological model that can accurately reproduce the vertical mixing processes, and the resulting effects on the nutrient transport and primary production. To fulfill this goal, the 1-D physical submodel needs to compute wind, thermal and tidal mixing explicitly and simultaneously. [6] Here a 1-D physical model [Mellor, 2001] that includes a turbulence closure model of level 2.5 type developed by Mellor and Yamada [1982] and improved by Mellor [2001] is used to provide realistic mixing coefficients at each of the fine vertical layers for both the physical and biological models. We introduced a tidal forcing term calculated with tidal current harmonic constants derived from the mooring acoustic Doppler current profiler (ADCP). The biological submodel is adapted from Eslinger et al. [2001], a 1-D model in Prince William Sound, Alaska, with nine compartments: two phytoplankton groups (diatom and flagellates), three zooplankton (micro, small and large zooplankton), three nutrients (nitrogen, ammonium, silicon) and detritus. This model has more compartments than the Eslinger and Iverson [2001] 1-D model for the Bering Sea, which has only four compartments: phytoplankton, zooplankton, nutrient and detritus. [7] The model was applied to the mooring site M2 (Figure 1). The following meteorological and ocean observations from the mooring are used in this study: wind speed, air temperature, fluorescence, current, ocean temperature and salinity at certain depth. Because mooring data are usually limited in both time and space, meteorological data sets, such as NCEP and ECMWF with global and daily or 6-hourly coverage, are commonly used as forcing in the long-period model runs and 3-D models. NCEP meteorological data is used here as model forcing for comparisons with mooring forced model results. A series of sensitivity studies were conducted to examine the phytoplankton bloom responses to the model parameters and processes. [8] In section 2, tidal current and meteorological data at M2 biophysical mooring site are introduced. The coupled physical-biological model equations and numerical experiment set up are presented in sections 2 to 4. Model results are discussed in section 5, followed by conclusions and discussions.

2. Physical Model 2.1. Physical Model Equations [9] The physical model includes the following onedimensional prognostic equations of horizontal velocity (U, V), potential temperature (T), salinity (S), and turbulence

2 of 20


C03002

Table 1. Magnitude of Tidal Harmonic Constants of Elevation and Depth-Averaged Current Component at the Mooring Site Elevation, cm u, cm/s v, cm/s

M2

K1

O1

N2

S2

Residual

45.42 19.65 14.06

24.96 13.87 1.91

17.02 8.15 0.82

14.01 4.05 2.54

0.76 0.75 0.58

2.88 0.07

components (q2, q2l) from a new version of the M-Y 2.5level turbulence model [Mellor, 2001], DU 1 @p @ @U ; fV þ ¼ KM Dt r0 @x @z @z

ð1Þ

DV 1 @p @ @V ; þ fU þ ¼ KM Dt r0 @y @z @z

ð2Þ

DT @ @T @R þ ¼ KH ; Dt @z @z @z

ð3Þ

DS @ @S ; ¼ KH Dt @z @z

ð4Þ

" # Dq2 @ @q2 @U 2 @V 2 2g @~r þ 2KM þ KH þ ¼ Kq Dt @z @z @z @z r0 @z 2e;

ð5Þ

( " # Dq2 l @ @q2 l @U 2 @V 2 þ ¼ Kq þ E1 l K M Dt @z @z @z @z ) g @~r ~: þ E3 K H leW r0 @z

ð6Þ

[10] Here f is the Coriolis parameter (s1); r0 is reference density; p is pressure; @R/@z is short wave radiation flux q3 ~ , W 1 + E2 (l/kL) where L1 jzj1 + divergence. e B1 l jD + zj1 is a ‘‘wall proximity’’ function and 0 < z < D.Vertical kinematic viscosity and vertical diffusivity are expressed as KM ¼ qlSM þ mmol

ð7aÞ

KH ¼ qlSH þ mmol ;

ð7bÞ

where mmol = 2.3 105 m2/s is the molecular diffusion coefficient, which has been adjusted to be slightly larger than the one used in POM (mmol = 2.0 105 m2/s) by comparing model results with mooring data. The coeffi-

C03002

cients, SM and SH, are functions of the Richardson number l 2 g @~r ) given by (GH = 2 q r0 @z SH ½1 ð3A2 B2 þ 18A1 A2 ÞGH ¼ A2 ½1 6A1 =B1

SM ½1 9A1 A2 GH SH 18A21 þ 9A1 A2 GH ¼ A1 ½1 3C1 6A1 =B1 :

ð8aÞ

ð8bÞ

[11] The factor, @r/@z, is the vertical density gradient minus the adiabatic lapse rate. KBqB is the vertical turbulence diffusivity which is practically set equal to KH in the simulation. (A1, B1, A2, B2, C1, E1, E2, E3) = (0.92, 16.6, 0.74, 10.1, 0.08, 1.8, 1.33, 1.0) are constants [see Mellor and Yamada, 1982]. Details of the numerical solutions of the model are given by Mellor [2001, 2003]. 2.2. Tidal Forcing [12] Since biological mechanisms are basically posed as purely vertical, the advantage of including tidal forcing in a 1-D physical model is to avoid some artificial mixing caused by parameterizations used by other 1-D models, such as Eslinger and Iverson [2001], Eslinger et al. [2001] and Merico et al. [2004]. Thus this can be an efficient and natural way to examine tidal mixing effects on an ecosystem. [13] ADCP measurements at the mooring site are available at depths between 10 m and 62 m every 2 m from 25 April to 18 June 2000. Tidal current harmonic analysis was conducted at each depth with five tidal constituents (M2, S2, N2, O1, K1) and residual current. Note that the mooring site name, also called M2, should be distinguished from the tidal constituent name M2. The tidal elevation harmonic constants can be derived from global geosat altimetry [Cartwright and Ray, 1990] at the mooring site. The magnitude of the tidal harmonic constants of elevation and the depth-averaged current component are listed in Table 1. The MB2 component is dominant and followed by K1, O1, N2 and S2, which is consistent with previous tidal studies in the region [Kowalik and Stabeno, 1999; Stabeno and van Meurs, 1999; Hermann et al., 2002]. [14] The vertical profiles of the tidal current constituents are shown in Figure 2 (the depth-averaged value was plotted at depth 0 m). The M2, N2 constituents and residual current had a peak at around 15 m, while the other tidal constituents were vertically homogenous. The semidiurnal tidal constituents M2 and N2 have thicker bottom boundary layers than the diurnal constituents because the inertial period (2p/f = 14.32 hr) at the mooring latitude is close to the semidiurnal tidal period. Thus semidiurnal tides contribute more to vertical mixing owing to their large vertical shear. The depth-averaged residual or long-term current (U0, V0) = (2.88cm/s, 0.07cm/s), larger than1 cm/s estimated by Coachman [1986], is still 1 order of magnitude smaller than total depth-averaged current of 25 to 45 cm/s. A small mean flow supports the assumption that it is reasonable to use a 1-D vertical model in the region. [15] With the depth-averaged tidal current harmonic constants obtained from the ADCP, tidal forcing (Tx, Ty) in a

3 of 20


C03002

KH ð@T =@z; @S=@zÞjz¼0 ¼ ðQlw þ Qsh þ Qlh Þ= Cp rw ;

ð E PÞSz¼0 ;

Figure 2. Vertical profiles of the five tidal constituents (M2, S2, N2, O1, K1) and residual current.

depth-averaged barotropic tide can be obtained from the equations (1) and (2), qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 @p DUT fVT þ Cd*UT * UT2 þ VT2 =H ¼ Dt r0 @x

ð9Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 @p DVT þ fUT þ Cd*VT * UT2 þ VT2 =H; ¼ Dt r0 @y

ð10Þ

Tx ¼

Ty ¼

where (UT, VT) are depth-averaged tidal current calculated from tidal current harmonic constants obtained from the ADCP, ðUT ; VT Þ ¼ ðU0 ; V0 Þ þ

X

Fn ðUMn ; VMn Þ cos½sn t þ Pn

n

ðUPn ; VPn Þ ;

ð11Þ

where Fn, sn, Pn, (UMn, VMn), (UPn, VPn) are the nodal modulation factor for amplitude, frequency, sum of the nodal and astronomical adjustment factors for phase; tidal current magnitude and phase lag (or tidal current harmonic constants) of the nth tidal current constituents, respectively. Cd is bottom drag coefficient typical in 2-D shallow water equations. Since the observed tidal current harmonic constants only cover 10 –62 m out of the total water depth H = 74 m, the value of Cd will be discussed later for Cd = 0 and Cd = 103, a typical value. 2.3. Surface Flux [16] Surface boundary includes wind stress, heat and salt flux, rw KM ð@U =@z; @V =@zÞjz¼0 ¼ tx ; ty ! ! ¼ rA CD U 10 U ðU10 U ; V10 V Þ ð12Þ

C03002

ð13Þ

where rw and rA are seawater and air density. Cp = 4190 J (Kg K)1 is specific heat of seawater. U10 and V10 are wind velocity at 10 m. CD is air-water drag coefficient. Qlw, Qsh, Qlh are net long wave radiation, sensible heat flux and latent heat flux, respectively. E and P are evaporation and precipitation. The sensible and latent heat flux, evaporation were calculated using bulk formulas by Gill [1982]. The long-wave and shortwave radiation are calculated following Rosati and Miyakoda [1988] and Parkinson and Washington [1979]. [17] The following meteorological variables are needed in the calculation of equations (13): wind speed, cloud cover, air temperature, precipitation rate, sea level pressure, and specific (or relative) humidity. All of these variables are available from NCEP reanalysis data at 6-hourly time interval provided by the NOAA-CIRES Climate Diagnostics Center, Boulder, Colorado, USA, from their Web site at http://www.cdc.noaa.gov/. A problem list is also updated on the Web with known deficiencies in the NCEP reanalysis data. forcing. However, at the M2 mooring site, only wind speed, air temperature and sea level pressure are available with good quality during 25 April to 20 September of 2000. An alternative solution for equations (13) is to directly prescribe sea surface temperature (SST) and sea surface salinity (SSS) from mooring data. The latter can provide a precise solution for the surface flux, but is usually limited by the time duration and spatial coverage of the observations. Thus both methods are used for the sensitivity studies.

3. Biological Model [18] The biological model was adapted from Eslinger et al. [2001], and the code was rewritten in a 3-D framework [Wang et al., 2003]. The model has nine compartments: two phytoplankton (diatom and flagellates: D and F), three zooplankton (small copepods, large copepods, and microzooplankton: ZS, ZL, ZP), three nutrients (nitrate+nitrite, ammonium, silicon: NO3, NH4, Si) and detritus (Det). The linkages between the nine compartments are illustrated in Figure 3. Note we have changed one large zooplankton type from Eslinger et al. [2001] to be microzooplankton, according to the fact that large zooplankton were rarely observed in the southeastern Bering Sea middle shelf. This biological submodel is coupled with the physical submodel described above. The vertical mixing coefficient KH is used for all the biological variables according to the conventional understanding that a planktonic organism moves passively with the water. The vertical advection term is used only for compartments that can sink by itself, such as for phytoplankton and detritus. The nutrients are neutral in water so no sinking term is used. Tidal mixing is explicitly calculated in the model instead of the tidal parameterization of Eslinger et al. [2001]. The biological model equations are expressed as follows: @D ¼ D GD RD RgD GDS ZS GDL ZL GDP ZP @t @ ð W D DÞ @ @D ; þ þ KH @z @z @z

4 of 20

ð14Þ

C03002


C03002

Figure 3. Diagram of linkages between the nine biological compartments. @F ¼ F GF RF RgF GFS ZS GFL ZL GFP ZP @t @ ðW F F Þ @ @F ; þ þ KH @z @z @z

ð15Þ

@ @ZS @ZS ; ¼ ZS AS GDS þ GFS 1 ExS M S þ KH @t @z @z ð16Þ

@ðW L ZLÞ @ZL ¼ ZL AL GDL þ GFL 1 ExL M L þ @t @z @ @ZL ; þ KH @z @z

@Det ¼ ZS 1 AS GDS þ GFS þ M S þ ZL @t

1 AL GDL þ GFL þ M L þ ZP

1 AP GDP þ GPM þ M P Det RgDet @ ðW Det DetÞ @ Det þ þ KH ; @z @z @z

where superscripts D and F denote diatoms and flagellates. Superscripts S, L, and P denote small copepods, large copepods, and micro-zooplankton. Terms G, R, Rg are phytoplankton growth rate, respiration rate, and mortality rate, respectively.

ð17Þ

@ @ZP @ZP ; ¼ ZP AP GDP þ GFP 1 ExP M P þ KH @t @z @z ð18Þ

@NO3 ¼ fNO3 D GD RD þ F GF RF þ CAtoN NH4 @t @ @NO3 þ ; ð19Þ KH @z @z

@NH4 ¼ ZS AS GDS þ GFS ExS þ ZL AL GDL þ GFL ExL þ ZP @t AP GDM þ GFM ExM þ D RgD þ F RgF þ Det

RgDet ð1 fNO3 Þ D GD RD þ F GF RF @ @NH4 ; ð20Þ CAtoN NH4 þ KH @z @z @ @Si @Si ; ¼ kSi D GD RD þ KH @t @z @z

ð21Þ

ð22Þ

GX ¼ mX0 erT min Nfrac ; Sifrac ; Ifrac ;

ð23aÞ

RX ¼ 0:05mX0 erT ;

ð23bÞ

RgX ¼ Rg0 erg T :

ð23cÞ

The superscript X denotes D or F. mX0 is maximum phytoplankton growth rate at 0C as discussed later. Nfrac, Sifrac, Ifrac are unitless ratios expressing nitrogen, silicon and light limitation as defined by Eslinger et al. [2001]. Sifrac is only used for diatoms. GXY denote the grazing rate of zooplankton Y on phytoplankton X, formulated as modified Ivlev-type grazing [Ivlev, 1945; Eslinger et al., 2001], GXY ¼ GXY 1 elXY ð X XY Þ :

[19] The sinking rate of phytoplankton and zooplankton (WX, WY), the remineralization rate of detritus (RgDet), and the ratio of phytoplankton growth due to nitrate uptake to that due to both nitrate and ammonium uptake (fNO3) follow Eslinger et al. [2001]. Table 2 provides a list of parameters used in the model, their units, and values. Further modifi-

5 of 20

C03002


C03002

Table 2. List of Parameter Values Symbol (Definition)

Unit 1

Rg0 (maximum phytoplankton mortality rate at 0C) r, rg (temperature coefficients) GDS, GDL, GDM, GFS, GFL, GFM (Ivlev-type grazing constants) lDS, lDL, lDM, lFS, lFL, lFM (Ivlev-type grazing constants) DS, DL, DM, FS, FL, FM (Ivlev-type grazing constants) AY (assimilation ratio of grazed phytoplankton) ExY (egestion ratio of grazed phytoplankton) MY (natural mortality rate of zooplankton) kSi (silicon to nitrogen mass ratio for diatoms)

h deg1 h1 mM1 mM N

d1

cations and improvements over Eslinger’s model are summarized below. [20] 1. The vertical mixing coefficient KH is used in each equation. Symbol for detritus has changed from Fec to Det here. [21] 2. A nitrification term is included in the nitrogen and ammonium equations. Ammonium is converted to nitrogen at a rate of CAtoN = 6.25 104 h1 (=1.5% d1). [22] 3. Silicon losses due to uptake are restored back to the bottom layer instantly, similar to Chai et al. [2002], who restored silicate below euphotic layer to its initial conditions. This is important to conserve the nutrients so that the model can run longer than one season. [23] 4. Maximum phytoplankton growth rate at 0C used F in equations (23a) and (23b), (mD 0 , m0 ) for diatom and 1 flagellate, were both 0.06 h in the work of Eslinger et al. [2001], and were 1.2 d1 and 0.65 d1, respectively (=0.05 h1, 0.027 h1) in the work of Merico et al. [2004]. Observations in the southeastern Bering Sea indicate that diatoms dominate spring bloom. Thus, in this study, we F used the same mD 0 to m0 ratio as used by Merico et al. [2004], but chose to increase the values slightly to 0.07 h1 and 0.038 h1 through model-data comparisons using different F values of mD 0 and m0.

4. Numerical Experiments [24] The 1-D model was applied to the mooring site M2 (Figure 1). The water depth H = 74 m. There are 37 vertical layers with 2 m per layer for both the physical and biological models. Model time step is 2 min. [25] The physical model was forced by tides, wind, shortwave radiation, surface heat and salt flux or prescribed SST and SSS. Different types of forcing were used to study the model sensitivity and the influence of the physical environment on the biological output. The standard run was designed to use observed forcing as much as possible, and thus was driven by tidal forcing, mooring wind, shortwave radiation, SST and SSS. The simulations run from 26 April to September 10 of 2000, when both the meteorological and oceanic mooring data were available and of good quality. The water temperature observations were available at 15 vertical layers (1, 6, 8, 12, 15, 20, 24, 27, 31, 34, 39, 44, 50, 56, 62 m). Salinity observations were available at seven vertical layers (1, 6, 12, 15, 24, 44, 62 m). Temperature and salinity at 1 m (or 6 m if data at 1 m were missing, which mainly occurred in salinity data) were used

Value 9.23 104 0.0633, 0.03 0.005, 0.003, 0.01, 0.013, 0.003, 0.01 1.50, 2.73, 2.73, 2.73, 2.73, 2.73 0.25, 0.50, 0.50, 0.25, 0.50, 0.50 0.7 0.4 0.06 1.104

as prescribed SST and SSS in the physical model. The alternative runs using NCEP meteorological data to calculate surface heat and salt flux were used in model sensitivity studies for comparison with the standard run. [26] The NCEP air temperature, wind speed, and calculated shortwave radiation using NCEP meteorological data are compared with the mooring data in Figure 4. The mean values and differences of the two data sets during the mooring period are presented in Table 3. The NCEP air temperature is 0.57C colder on average, wind speed is 0.99m/s stronger than that of the mooring, and the calculated shortwave radiation is very close to the observation (5% difference). [27] The standard run and eight sensitivity cases were designed to study the influences of the forcing terms on both the physical environments and biological model results. The forcing for each case is described below and listed in detail in Table 4. [28] 1. Case 1 has no tidal forcing in the standard run. [29] 2. Case 2 has no wind forcing in the standard run. [30] 3. Case 3 has constant T and S in the standard run. Thus no stratification will be formed. [31] 4. Case 4 has NCEP wind replacing mooring wind from the standard run. [32] 5. Case 5 has calculated shortwave radiation replacing mooring shortwave radiation from Case 4. [33] 6. Case 6 uses heat flux to force model instead of prescribed SST from case 5. [34] 7. Case 7 uses salt flux to force model instead of prescribed SSS from Case 6. Thus no mooring data are used for forcing in this case. The case runs with the same time duration as the standard run. [35] 8. Case 8 is the same as Case 7, but runs for the full year of 2000. [36] Initial velocity was set to zero. Initial temperature and salinity were vertically homogeneous based on the mooring data: T = 0.43C, and S = 31.97, if model starts from April 26; T = 1.46C, and S = 31.97, if model started from January 1. [37] Initial conditions for the biological model were also vertically homogeneous based on historical measurements and National Oceanographic Data Center (NODC)’s World Ocean Atlas (WOA) 2001 climatology nitrate data: D = 0.2 mmol N m3, F = 0.05 mmol N m3, ZS = 0.3 mmol N m3, ZL = 0.01mmol Nm3, ZM = 0.01 mmol N m3, NO3 = 10 mmol N m3, NH4 = 0 mmol N m3, Si = 40 mmol Si m3, Det = 0 mmol N m3. The N-based unit can be

6 of 20

C03002


Figure 4. Comparison of NCEP and mooring forcing data: (a) air temperature, (b) wind speed, and (c) daily mean short wave radiation. converted to the chl-based unit assuming a carbon:chlorophyll mass ratio of 40 g C (g chl)1 [Eslinger and Iverson, 2001] and a C:N molar ratio of 6.625 [Redfield et al., 1963].

5. Results and Discussion 5.1. Physical Model Validation 5.1.1. Current [38] The depth-averaged tidal currents reconstructed with the harmonic constants and simulated with the 1-D model compare well with the ADCP data (Figure 5), which indicates that the analyzed tidal constituents are the major ones in the region. The model started from day 115 (26 April), and the figure represents days 120 to 129. Since the semidiurnal tidal current is dominant, the flow is adjusted to observations quickly in 5 days. Here only 9 days of data are shown for a clear view of daily variations of the current and the rest of the model days also were compared similar to the 9 days (not shown). The bottom drag coefficient Cd (used in equations (9) and (10)) has a value of 0.001 in typical 2-D tidal models. The results using Cd = 0 and Cd = 103 show only negligible differences in

C03002

the simulated thermocline structure and current speed (not shown). Thus the model is not sensitive to Cd, and Cd = 0 was used for all the simulations. 5.1.2. Temperature, Salinity, and MLD [39] The time series of the simulated and observed temperature and MLD are shown in Figures 6 and 7. Because salinity contains lots of missing data, especially in the middle layers, the simulated and observed salinity were only compared in the surface layer (6 m) and the bottom layer (44 m) in Figure 8. Modeled and observed temperature was compared in the surface layer (1 m) and bottom layer (62 m) in Figure 9. The MLD for both the observations and simulations is defined as the depth at which the temperature becomes 0.2C less than the SST, following Eslinger and Iverson [2001], Kuhn and Radach [1997], and Fujii et al. [2002]. Since salinity varied in a narrow range of 31.8 to 32 psu in the both the observations and the standard run (Figure 8), the main structure of thermocline was determined by temperature. [40] The formation of the thermohaline in early May, and the progression of the upper and bottom mixed layers, were well reproduced in the standard run (Figure 6). In contrast, the bottom mixed layer was not developed in case 1 (no tidal forcing) and the upper mixed layer was very shallow or almost nonexisting in case 2 (no wind forcing). The depth of the 4C contour in Figure 6 at the end of the model run were 29 m, 34 m, and 19 m for the standard run, case 1 and case 2, respectively, which indicates that wind mixing deepens MLD and tidal mixing deepens the bottom boundary layer. The tidal and wind mixing effects can be easily identified by comparing the vertical mixing coefficient in Figures 10a– 10c. The vertical mixing coefficient in case 8 (Figure 10d) illustrated the annual cycle of the mixed layer forming in the spring, deepening in the summer to fall and disappearing in the winter. [41] Case 4 with NCEP wind forcing had a deeper MLD around day 160 and 220 related to stronger NCEP wind speed than the mooring (Figures 4b and 7). Therefore the bottom temperature in case 4 progressed faster than that of the standard run (Figures 6b and 6e). A cooler air temperature of NCEP compared to the mooring (Table 3) also contributed to deeper MLD, as seen by comparing Figures 7d and 7f. [42] The simulated surface and bottom salinity and temperature of the standard run, case 7 and case 8 are shown in Figures 8 and 9. In the standard run, SST and SSS were prescribed with observations. The simulated bottom temperature agrees very well with the observations, but the bottom salinity missed the fluctuations in the observations. The rising salinity above 32 psu between days 200 and 220 (Figure 8) can be caused by advection due to shelf break water intrusion, because SSS had been under the initial

Table 3. Mean Values and Differences of Air Temperature, Wind Speed, and Shortwave Radiation From NCEP and Mooring Data During the Mooring Period

Air temperature, C Wind speed, m/s Shortwave radiation, W/m2

7 of 20

NCEP

Mooring

Difference (NCEP Mooring)

6.56 6.35 154

7.13 5.36 146

0.57 0.99 8


C03002

C03002

Table 4. Numerical Experiments Design for the Standard Run and Sensitivity Studiesa Tidal forcing Stratification Mooring wind NCEP wind Mooring SW Calculated SW Mooring SST Heat flux Mooring SSS Salt flux Model time 26 April to 10 Sept Model time 1 Jan 1 to 31 Dec

Standard Run

Case 1

Case 2

Case 3

Case 4

Case 5

Case 6

Case 7

Case 8

Y Y Y

Y Y

Y Y

Y

Y Y

Y Y

Y Y

Y Y

Y Y

Y

Y

Y

Y Y

Y

Y

Y

Y

Y

Y

Y

Y Y

Y

Y

Y

Y

Y

Y

Y

Y

Y

Y Y

Y

Y

Y

Y

Y

Y

Y

Y Y

Y

Y

Y Y

Y

Y

a

Y denotes yes. Blank entries denote ‘‘no’’ or ‘‘not applicable.’’

31.97 psu and declining with time. Thus the fluctuations of the bottom salinity may be caused by this cross-shelf transport from the inner domain with lower salinity between days 130 and 180 and the outer domain with higher salinity between days 200 and 220. The simulated SST in cases 7 and 8 followed the observed fluctuations very well and compared well to observations before day 150. But afterwards the simulated SST was lower than the observations, with the largest discrepancy of 2C, which was mainly caused by lower NCEP air temperature (Table 3). The simulated bottom temperature had only a slight discrepancy with observations. The simulated SSS followed with the mean observed SSS trend well and captured the fluctuations between days 200 and 260, but missed the fluctuations between days 130 and 160, when advections occurred. 5.2. Biological Model Validation [43] Simulated total phytoplankton biomass for all cases was compared with fluorescence measurements at 12 m, 24 m, and 44 m, respectively, as shown in Figure 11 (12-m data published in work by Stabeno et al. [2001]; 24-m and 44-m data are unpublished). [44] Fluorescence serves as a proxy for phytoplankton biomass and can be converted to units of mg Chl mP3P using an appropriate conversion factor, which for this data set has not yet been determined [Stabeno et al., 2001; Hunt and Stabeno, 2002]. For comparison purposes, the fluorom-

eter data were multiplied to bring the data to a similar magnitude as the total phytoplankton of the standard run. [45] The onset, magnitude and duration of the spring phytoplankton blooms observed in the fluorescence data were well reproduced in the model results at each layer of the standard run (Figure 11). The timing differences of the simulated and observed bloom maximum were 1 day at 12 m and 24 m, and 2 days at 44 m. The simulated post blooms after day 200 show good agreement in timing compared with the fluorescence data. The spring phytoplankton bloom was dominated by diatoms, and the post blooms by flagellates, which agree with previous studies [Sukhanova et al., 1999; Merico et al., 2004]. The peak phytoplankton concentration reached 8 mmol N m3, or 16 mg Chl m3, the same order of magnitude as observed in PROBES program in 1980 and 1981 [Eslinger and Iverson, 2001] in the southeastern Bering Sea middle shelf. [46] The annual cycle of the simulated nitrate in case 8 was similar to that of the WOA2001 nitrate climatology (Figure 12) from May to September, but different in winter months (January to March), since late-autumn to winter cross-shelf mixing and ice-associated activities were not included in this model. 5.3. Sensitivity Analysis [47] The important characteristics for sensitivity simulations are presented in Table 5. MMLD is defined as the

Figure 5. Comparison of depth-averaged velocity from the ADCP-observed, simulated in the standard run, and reconstructed using tidal current harmonic constants: (a) u and (b) v. 8 of 20

C03002


C03002

Figure 6. Temperature against time: (a) mooring, (b) standard run, (c) case 1: no tide, (d) case 2: no wind, (e) case 4: NCEP wind, and (f) case 7: all forcing from NCEP.

mean MLD over the model time period Tm (26 April to 10 September); MPC is the mean phytoplankton concentration over depth H and time Tm; NPP is the net primary production; and GPP the gross primary production. These terms are mathematically defined below as the integration

over the time period Tm with time t, and in the entire water column for z,

9 of 20

1 MMLD ¼ Tm

Z

ð MLDÞdt;

ð24Þ

C03002


Figure 7

10 of 20

C03002

C03002


MPC ¼

1 HTm

Z

NPP ¼ NPPD þ NPPF ; NPPD ¼ Z NPPF ¼ F GF RF dzdt;

GPP ¼ GPPD þ GPPF ; GPPD ¼ Z GPPF ¼ FGF dzdt

Figure 8. Simulated and observed salinity at 6 m and 44 m: (a) standard run, (b) case 7, and (c) case 8.

C03002

ð D þ F Þdzdt; Z

ð25Þ

D GD RD dzdt; ð26Þ

Z

DGD dzdt; ð27Þ

where superscript D and F denote diatom and flagellate. The ratio of each case to the standard run for MPC, NPP, and GPP, and the ratios of NPPD/NPPF and GPPD/GPPF are also included in Table 5 for the convenience of discussion. [48] Annual primary production in the Bering Sea middle shelf has been estimated from observations using numerous methods (e.g., nitrate disappearance, carbon uptake studies). The estimates of annual primary production on the Bering Sea middle shelf (Table 6) have a wide range of 60 to 180 g C m2, or 750 to 2250 mmol N m2. Here a conversion factor of 1 g C m2 = 12.5 mmol N m2 is used according to Redfield et al. [1963]. Most of the NPP and GPP from the standard run and cases 1 – 8 (Table 5) were within this range. Those estimates are best compared to the NPP from the model, indicating that the modeled productivity is on the low side. However, the low NPP from the model is probably due to initialization of the model with a nitrate concentration of 10 mM, which is at the low end of the 10- to 15-mM range displayed for the middle shelf during April in Figure 7 of Whitledge and Luchin [1999]. 5.3.1. Effects of Tidal, Wind Mixing, and Thermal Stratification [49] Comparison of the no-tide, no-wind, and no stratification cases (cases 1, 2, and 3) with the standard case revealed distinct effects of tidal mixing, wind stirring, thermal stratification and their interactions on the timing and magnitude of the phytoplankton bloom and the gross and net primary production (Figure 11 and Table 5). [50] The absence of tidal mixing (case 1) had little influence on the onset of the spring phytoplankton bloom at 12 m in the model, since the spring bloom was triggered by the easing of the wind and development of stratification. However, the bloom declined faster, and post blooms were weaker than that of the standard run. The onset of the bloom was shifted 1 week earlier at 24 m, and 5 days later at 44 m. MPC, GPP, and NPP decreased 20% from that of the standard run, because nutrient replenishment into the euphotic zone was diminished or lessened without tidal mixing.

Figure 7. Comparison of the simulated and observed mixed-layer depth: (a) standard run, (b) case 1: no tide, (c) case 2: no wind, (d) case 4: NCEP wind, (e) case 5: NCEP wind and shortwave radiation, (f) case 6: NCEP wind, shortwave radiation, and surface heat flux, (g) case 7: short run with all forcing from NCEP, and (g) case 7: 1-year run with all forcing from NCEP. 11 of 20

C03002


Figure 9. Simulated and observed temperature at 1 m and 62 m: (a) standard run, (b) case 7, and (c) case 8.

[51] Without wind mixing (case 2), the timing, duration and magnitude of the phytoplankton bloom changed significantly from the standard run. A prolonged phytoplankton bloom occurred at 12 m in the absence of wind mixing.

C03002

Note that since the MLD was very shallow without wind mixing, the 12 m depth was below the mixed layer in the spring and early summer, while at the thermocline layer in late summer. Although GPP and NPP decreased by 4% and 13%, respectively, MPC increased 20% (Table 5), probably owing to more phytoplankton being mixed to deep water where it was subjected to less zooplankton grazing. The simulated maximum phytoplankton biomass at 24 m and 44 m were both 2.6 mmol N m3P, different from the pattern in the standard run with 5 mmol N m3 at 24 m, and 1.85 mmol N m3 at 44 m. [52] With no stratification (case 3), the water column was vertically well mixed, and the phytoplankton blooms were similar in timing, duration and magnitude throughout the water column for a longer period than in the standard run. The simulated phytoplankton concentrations were almost identical at 12 m, 24 m and 44 m because of strong vertical mixing. Thus the concentration at 12 m was much lower, and at 44 m was higher than in the standard run. GPP increased by 8% in this case, probably because all the nutrients in the water column were available for photosynthesis. MPC was 45% more than that of the standard run owing to less zooplankton grazing in the bottom. However, NPP was 19% less, because more phytoplankton was mixed below the euphotic zone, where loss of phytoplankton biomass by respiration exceeds production by photosynthesis. 5.3.2. Effects of NCEP Versus Mooring Meteorological Forcing [53] By replacing the mooring wind, shortwave radiation, and surface heat flux with NCEP data (cases 4 to 6), the phytoplankton bloom was similar to the standard run (Figure 11), so, case 5 was not shown in the figure. By using salt flux in cases 7 and 8, the magnitude of the first phytoplankton bloom was reduced and the second peak was higher and lasted longer, although the post blooms were similar to the standard run. The simulated SSS between days 130 and 160 were smoother than the observations (Figure 8). This reduced MLD fluctuations and nutrient entrainment into the surface layer, resulting in nutrient limitation that changed the magnitude and timing of the bloom. [54] The results of a 1-year run with all forcing using NCEP data (case 8) were similar in timing and magnitude to the phytoplankton bloom of case 7, but the peak bottom phytoplankton concentration at 44m was slightly higher (Figure 11). MMLD was similar but the NPP and GPP were higher. The GPP and NPP for cases 7 and 8 were different because the initial physical-biological conditions for case 7 were slightly different from that of case 8 on 26 April. [55] The MMLD in the standard run and cases 4 to 7 (Table 5) displayed differences due to the bias in NCEP forcing (Table 3): stronger mean winds and lower air temperatures contributed to a deeper MMLD, while higher shortwave radiation (case 5) contributed to a shallower MMLD. MPC, GPP and NPP increased from cases 4 to 7 (Table 5), and were higher than that of the standard run. The increase from case 4 to 5 can be explained by the increased light for photosynthesis. The increase from the standard run to cases 4, 6 and 7 was related to the increase of the MMLD.

12 of 20


C03002

C03002

Figure 10. Simulated vertical mixing coefficient KH in (a) standard run, (b) case 1, (c) case 2, and (d) case 8.

5.3.3. Linkage of MLD to Primary Production and Phytoplankton Groups [56] To investigate the relationship between primary production and MLD, we designed an experiment using the standard run, but with a fixed vertical mixing coefficient profile for the biological model only, KH ¼

0:1 m2 =s; 2 105 m2 =s;

z < MLD; or z > MLD þ 4 MLD < z < MLD þ 4:

[57] MLD represents a fixed mixed-layer depth for each simulation. The two kinds of the mixing coefficients, 0.1 and 2 105 m2/s, denote turbulent mixing and molecular mixing, respectively. There is a 4-m-thick pycnocline between the surface and bottom mixed layers. We have run cases with MLD varying from 5 m to 60 m. GPP and NPP versus MLD are shown in Figure 13a. There exists a maximum GPP at MLD = 24 m, while NPP did not

show a similar trend. This idealized MLD hypothesis can explain the changes of GPP with MMLD from standard run to cases 4, 6 and 7. However, the NPP may be affected more by the history of MLD variation than by the MMLD. [58] Figure 13b reveals that the dominant phytoplankton group in both GPP and NPP is a function of MLD. Flagellates dominate when MLD < 15m, while diatoms dominate when MLD > 15 m. The ratio of NPPD/NPPF and GPPD/GPPF for the standard run and cases 2 to 7 (Table 5) did show a correspondent increase with MLD, with the smallest ratio in case 2 (no-wind, shallowest MLD) and the highest in case 3 (no stratification, MLD is infinite or equal to the water depth). Around MMLD = 14 m, the NPPD/ NPPF ratio switched from flagellate dominant to diatom dominant; nevertheless for GPP, diatoms were dominant in all cases. Thus the ratio of GPPD/GPPF may be affected more by the history of MLD variation than by the MMLD.

13 of 20

C03002


Figure 11. Time series of fluorometer and simulated total phytoplankton, diatom, flagellate, and total zooplankton concentrations at 12 m, 24 m, and 44 m in the standard run and cases 1 to 8, except case 5. 14 of 20

Figure 11. (continued)

C03002

C03002




15 of 20

C03002

C03002




16 of 20

C03002

C03002




17 of 20

C03002


C03002

Figure 12. Annual cycle of the nitrate from (a) WOA2001 and (b) case 8.

Case 1 is not explainable using the above MLD theory because the lack of tidal mixing to bring bottom nutrient up had significantly reduced its primary production.

6. Conclusions and Discussion [59] A coupled physical-biological model was developed for the southeastern Bering Sea middle shelf. A

C03002

tidal forcing term was introduced into the 1-D physical model to achieve a realistic presentation of tidal mixing instead of parameterizations of tidal mixing used in the previous 1-D biological models [Eslinger et al., 2001]. By using a turbulence closure model of level 2.5 type [Mellor, 2001], tidal mixing and its interactions with wind stirring, and thermal stratification were realistically reproduced that provided an important tool for investigating the implications of these mixing terms and their interactions on the biological processes. We improved a 1-D biological model [Eslinger et al., 2001] by introducing vertical mixing coefficient KH from the turbulence closure model, adding nutrient regeneration sources to the equations, and setting the maximum phytoplankton growth F rate at 0C, (mD 0 , m0 ), to suit the ecosystem of the southeastern Bering Sea middle shelf. [60] The 1-D coupled physical-biological model successfully reproduced the observed currents, thermocline, and MLD, and the onset, magnitude and duration of the spring phytoplankton blooms at 12 m, 24 m, and 44 m. The spring phytoplankton bloom was dominated by diatoms, and the post blooms by flagellates, which agree with previous studies in the region. The peak phytoplankton concentration reached 8 mmol N m3, or 16 mg Chl m3, the same order as was observed in the PROBES program in 1980 and 1981 [Eslinger and Iverson, 2001]. The NPP and GPP of the standard run and cases 1 – 8 were within the range of 60 to 200 g C m2(=750 to 2500 mmol N m2) estimated in previous studies. Sensitivity studies revealed distinct effects of tidal mixing, wind stirring, thermal stratification and their interactions on the timing and magnitude of the phytoplankton bloom and the gross and net primary production. Without tidal mixing, MPC, GPP, and NPP all decreased 1920% from that of the standard run. Without wind mixing, the timing, duration and magnitude of the phytoplankton bloom changed significantly from the standard run. Although GPP and NPP decreased 4% and 13%, respectively, MPC increased 20% due to reduced zooplankton grazing in the bottom mixed layer. With no stratification, the phytoplankton bloom can be seen throughout the water column in a longer period than in the standard run. GPP and MPC were 8% and 45% higher, respectively, but NPP was 19% less. Wind stirring and thermal stratification were critical to the onset and duration of the spring phytoplankton bloom, and tidal mixing was critical to the vertical nutrient transport to sustain the bloom and the annual primary production. [61] The NCEP meteorological data has a bias compared with the mooring data. The NCEP air temperature is 0.57C

Table 5. Important Characteristics for Simulations From 26 April to 10 September 2000a Standard Run MMLD, m MPC, m mol N m3 NPP, m mol N m2 NPPD/NPPF GPP, m mol N m2 GPPD/GPPF

12.14 1.028 100% 792 100% 0.82 1741 100% 1.86

Case 1

Case 2

Case 3

Case 4

Case 5

Case 6

Case 7

Case 8

12.31 0.829 81% 637 80% 0.65 1387 80% 1.57

5.34 1.235 120% 688 87% 0.53 1692 96% 1.14

None 1.493 145% 640 81% 4.5 1876 108% 6.76

14.01 1.071 104% 808 102% 1.03 1825 105% 2.26

13.80 1.102 107% 841 106% 1.01 1885 108% 2.21

16.78 1.143 111% 869 110% 1.17 1958 112% 2.53

19.69 1.229 119% 875 110% 1.41 2057 118% 3.03

19.11 1.365 133% 995 126% 0.79 2136 123% 1.55

a MMLD, MPC, NPP, and GPP are the important characteristics. Values in the second line for MPC, NPP, and GPP are the ratios of the results for each case to the results for the standard run. NPPD/NPPF and GPPD/GPPF represent the diatom to flagellate ratio.

18 of 20


C03002

C03002

Table 6. Estimates of the Annual Primary Production on the Bering Sea Middle Shelf and Related References Annual Primary Production, g C m2 150 – 175 60 – 180

75 – 125 63 139

References Sakshaug [2004], based on the work of Sambrotto et al. [1986], Walsh and McRoy [1986], Whitledge et al. [1986], Walsh and Dieterle [1994] Springer et al. [1996], based on the work of Codispoti et al. [1982], Sambrotto et al. [1984], Walsh and McRoy [1986], Whitledge et al. [1986], Hansell et al. [1993], Walsh and Dieterle [1994] Hansell et al. [1993] NPP by the standard run of this model study GPP by the standard run of this model study

colder, wind speed is 0.99m/s stronger than that of the mooring, and the calculated shortwave radiation is 5% higher than the observations. Sensitivity studies demonstrated that these biases can affect the MLD and phytoplankton bloom. Linkage of MLD to primary production and dominant phytoplankton species were found. If a constant MLD was used in the model, there exists a maximum GPP when MLD = 24 m. The dominate species of both the gross and net primary production were flagellates when MLD < 15 m, and diatoms when MLD > 15 m. This idealized hypothesis can explain the changes of GPP with MMLD from standard run to cases 4, 6 and 7. The NPPD/NPPF ratio switched around MMLD = 14 m from flagellate dominant to diatom dominant, but for GPP, diatom was dominant in all cases. Thus the ratio of GPPD/GPPF may be affected more by the history of MLD variation than by the MMLD. Field measurements are needed to validate the above hypothesis from the model studies. [62] The above studies suggest that errors in the physical environment can be magnified in the biological production mainly through the vertical mixing processes. This explains why the Bering Sea ecosystems are sensitive indicators of the effects of climate changes in the subarctic region [Stabeno et al., 2001; Overland and Stabeno, 2004]. Thus an accurate representation of all the mixing mechanisms (tide, wind, thermodynamics, etc.) in the physical model is critical to the success of the biological model, the same conclusion that many other coupled physical-biological modelers [i.e., Eslinger and Iverson, 2001; Hood et al., 2003] have also come to. Lack of an explicit mixed layer model can cause deficiencies in biological modeling [Kawamiya and Oschlies, 2003]. Because of lack of tidal forcing in the previous 1-D models, tidal mixing process has not been realistically included to interact with wind mixing and thermal mixing/stratification, i.e., models with two layers [Eslinger and Iverson, 2001; Merico et al., 2004] in shallow waters, models with multilayers [Hood et al., 2003; McCreary et al., 2001] in deep waters with both seasonal and permanent thermoclines, models using turbulence closure model [Kuhn and Radach, 1997; Denman and Pena, 1999]. We recommend using turbulence closure model because it can produce continuously varying profiles of mixing coefficients rather than step-like profiles in layered MLD models. This study provided a way to include a tidal forcing term in a vertically resolved 1-D model and demonstrated the importance of tidal mixing and its interaction with wind mixing and thermal mixing/

stratification on modeling MLD and biological processes in the Bering Sea shelf. Thus it is highly recommended to include tidal mixing explicitly in the biological modeling in the Bering Sea shelf and similar regions with strong tidal mixing.

Figure 13. Mixed-layer depth (MLD) against primary production: (a) GPP-D, GPP-F, NPP-D, and NPP-F and (b) GPP, NPP.

19 of 20

C03002


[63] Acknowledgments. We are grateful to Phyllis Stabeno and Sigrid Sallo, PMEL/NOAA, Seattle, Washington, for providing the M2 buoy and subsurface mooring data to drive and validate the 1-D model, and we thank NPRB and PMEL for supporting the collection of data at site 2. We also thank D. Eslinger for providing the ecosystem model code and helping on the model development. This study was supported by International Arctic Research Center (IARC), University of Alaska Fairbanks through the JAMSTEC-IARC Research Agreement.

References Cartwright, D. R., and R. D. Ray (1990), Oceanic tides from Geosat altimetry, J. Geophys. Res., 95, 3069 – 3090. Chai, F., R. C. Dugdale, T. H. Peng, F. P. Wilkerson, and R. T. Barber (2002), One-dimensional ecosystem model of the equatorial Pacific upwelling system: Part I. Model development and silicon and nitrogen cycle, Deep Sea Res., Part II, 49, 2713 – 2745. Coachman, L. K. (1986), Circulation, water masses, and fluxes on the southeastern Bering Sea shelf, Cont. Shelf Res., 5, 23 – 108. Codispoti, L. A., G. E. Friederich, R. L. Iverson, and D. W. Hood (1982), Temporal changes in the S.E. Bering Sea’s inorganic carbon system during Spring 1980, Nature, 296, 242 – 245. Denman, K. L., and M. A. Pena (1999), A coupled 1-D biological/physical model of the northeast subarctic Pacific Ocean with iron limitation, Deep Sea Res., Part II, 46, 2877 – 2908. Eslinger, D. L., and R. L. Iverson (2001), The effects of convective and wind-driven mixing on spring phytoplankton dynamics in the southeastern Bering Sea middle shelf domain, Cont. Shelf Res., 21, 627 – 650. Eslinger, D. L., R. T. Cooney, C. P. McRoy, A. Ward, T. Kline, E. P. Simpson, J. Wang, and J. R. Allen (2001), Plankton dynamics: Observed and modeled responses to physical forcing in Prince William Sound, Alaska, Fish. Oceanogr., 10, suppl. 1, 81 – 96. Francis, R. C., S. R. Hare, A. B. Hollowed, and W. S. Wooster (1998), Effects of interdecadal climate variability on the oceanic ecosystems of the NE Pacific, Fish. Oceanogr., 7, 11 – 21. Fujii, M., Y. Nojiri, Y. Yamanaka, and M. J. Kishi (2002), A one-dimensional ecosystem applied to time-series Station KNOT, Deep Sea Res., Part II, 49, 5441 – 5461. Gill, A. E. (1982), Atmosphere-Ocean Dynamics, 662 pp., Elsevier, New York. Hansell, D. A., T. E. Goering, and J. J. Walsh (1993), Patterns of nitrate utilization and new production over the Bering-Chukchi Shelf, Cont. Shelf Res., 13, 601 – 627. Hermann, A. J., P. J. Stabeno, D. B. Haidvogel, and D. L. Musgrave (2002), A regional tidal/subtidal circulation model of the southeastern Bering Sea: Development sensitivity analyses and hindcasting, Deep Sea Res., Part II, 49, 5945 – 5967. Hood, R. R., K. E. Kohler, J. P. McCreary, and S. L. Smith (2003), A fourdimensional validation of a coupled physical-biological model of the Arabian Sea, Deep Sea Res., Part II, 50, 2917 – 2945. Hunt, G. L., and P. J. Stabeno (2002), Climate change and the control of energy flow in the southeastern Bering Sea, Prog. Oceanogr., 55, 5 – 22. Ivlev, V. S. (1945), The biological productive of waters, Usp. Sovrem. Biol., 19, 98 – 102. Kawamiya, M., and A. Oschlies (2003), An eddy-permitting coupled ecosystem-circulation model of the Arabian Sea, J. Mar. Syst., 38, 221 – 257. Kowalik, Z., and P. Stabeno (1999), Trapped motion around the Pribilof Island in the Bering Sea, J. Geophys. Res., 104(C11), 25,667 – 25,684. Kuhn, W., and G. Radach (1997), A one-dimensional physical-biological model study of the pelagic nitrogen cycling during the spring bloom in the northern North Sea (FLEX’ 76), J. Mar. Res., 55, 687 – 734. McCreary, J. P., K. E. Kohler, R. R. Hood, S. Smith, J. Kindel, A. S. Fischer, and R. A. Weller (2001), Influences of diurnal and intraseasonal forcing on mixed-layer and biological variability in the central Arabian Sea, J. Geophys. Res., 106(C4), 7139 – 7155. Mellor, G. L. (2001), One-dimensional, ocean surface layer modeling: A problem and a solution, J. Phys. Oceanogr., 31, 790 – 809. Mellor, G. L. (2003), User’s guide for a three-dimensional primitive equation, numerical ocean model, report, Program in Atmos. and Ocean. Sci., Princeton Univ., Princeton, N. J. Mellor, G. L., and T. Yamada (1982), Development of a turbulent closure model for geophysical fluid problems, Rev. Geophys., 20, 851 – 875. Merico, A., T. Tyrrel, E. J. Lessard, T. Oguz, P. J. Stabeno, S. I. Zeeman, and T. E. Whitledge (2004), Modeling phytoplankton succession on the Bering Sea shelf: Role of climate influences and trophic interactions in generating Emiliania huxleyi blooms 1997 – 2000, Deep Sea Res., Part I, 51, 1803 – 1826. Napp, J. M., and G. L. Hunt Jr. (2001), Anomalous conditions in the southeastern Bering Sea, 1997: Linkages among climate, weather, ocean, and biology, Fish. Oceanogr., 10, 61 – 68.

C03002

Niebauer, H. J., V. A. Alexander, and S. Hendricks (1990), Physical and biological oceanographic interaction in the spring bloom at the Bering Sea marginal ice edge zone, J. Geophys. Res., 95(C12), 22,229 – 22,242. Niebauer, H. J., V. Alexander, and S. M. Henrichs (1995), A time-series study of the spring bloom at the Bering Sea ice-edge: I. Physical processes, chlorophyll, and nutrient chemistry, Cont. Shelf Res., 15, 1859 – 1877. Overland, J. E., and P. J. Stabeno (2004), Is the climate of the Bering Sea warming and impacting the ecosystem?, Eos Trans. AGU, 85(33), 309 – 310, 312. Overland, J. E., N. A. Bond, and J. Miletta (2001), North Pacific atmospheric and SST anomalies in 1997: Links to ENSO?, Fish. Oceanogr., 10, 69 – 80. Parkinson, C. L., and W. M. Washington (1979), A large-scale numerical model of sea ice, J. Geophys. Res., 84(C1), 311 – 337. Pollard, R. T., P. B. Rhines, and R. O. R. Y. Thompson (1973), The deepening of the wind-mixed layer, Geophys. Fluid Dyn., 4, 381 – 404. Redfield, A. C., B. H. Ketchum, and F. A. Richards (1963), The influence of organisms on the composition of seawater, in The Sea, vol. 2, edited by M. N. Hill, pp. 26 – 79, John Wiley, Hoboken, N. J. Rosati, A., and K. Miyakoda (1988), A general-circulation model for upperocean simulation, J. Phys. Oceanogr., 18, 1601 – 1626. Sakshaug, E. (2004), Primary and secondary production in the Arctic seas, in The Organic Carbon Cycle in the Arctic Ocean, edited by R. Stein and R. W. Macdonald, pp. 57 – 81, Springer, New York. Sambrotto, R. N., J. J. Goering, and C. P. McRoy (1984), Large yearly production ofphytoplankton in the western Bering Strait, Science, 225, 1147 – 1150. Sambrotto, R. N., H. J. Niebauer, J. J. Goering, and R. L. Iverson (1986), Relationships among vertical mixing, nitrate uptake and phytoplankton growth during the spring bloom in southeastern Bering Sea, Cont. Shelf Res., 5, 161 – 198. Schumacher, J. D., and V. Alexander (1999), Variability and influence of the physical environment to the ecosystem in the Bering Sea, in The Bering Sea: A Summary of Physical, Chemical and Biological Characteristics and a Synopsis of Research, edited by T. R. Loughlin and K. Ohtani, pp. 147 – 160, Alaska Sea Grant Press, Fairbanks. Springer, A. M., P. C. McRoy, and M. V. Flint (1996), The Bering Sea Green Belt: Shelf-edge processes and ecosystem production, Fish. Oceanogr., 5, 205 – 223. Stabeno, P. J., and P. van Meurs (1999), Evidence of episodic on-shelf flow in the southeastern Bering, J. Geophys. Res., 104(C12), 29,715 – 29,720. Stabeno, P. J., J. D. Schumacher, R. F. Davis, and J. M. Napp (1995), Under-ice observations of water column temperature, salinity and spring phytoplankton dynamics: Eastern Bering Sea shelf, J. Mar. Res., 56, 239 – 255. Stabeno, P. J., A. Bond, B. Kachel, A. Salo, and J. D. Schumacher (2001), On the temporal variability of the physical environment over the southeastern Bering Sea, Fish. Oceanogr., 10, 81 – 98. Sukhanova, I. N., H. J. Semina, and M. V. Venttsel (1999), Spatial distribution and temporal variability of phytoplankton in the Bering Sea, in Dynamics of the Bering Sea, pp. 453 – 484, Alaska Sea Grant Press, Fairbanks. Thompson, R. O. R. Y. (1976), Climatological numerical models of the surface mixed layer of the ocean, J. Phys. Oceanogr., 6, 496 – 503. Walsh, J. J., and D. A. Dieterle (1994), CO2 cycling in the coastal ocean: I. A numerical simulation analysis of the southeastern Bering Sea with applications to the Chukchi Sea and the northern Gulf of Mexico, Prog. Oceanogr., 34, 335 – 392. Walsh, J. J., and C. P. McRoy (1986), Ecosystem analysis in the southeastern Bering Sea, Cont. Shelf Res., 5, 259 – 288. Wang, J., C. Deal, Z. Wan, M. Jin, N. Tanaka, and M. Ikeda (2003), User’s guide for a physical-ecosystem model (PhEcoM) in the subpolar and polar oceans, version 1, Tech. Rep. 03-01, Int. Arctic Res. Cent., Fairbanks, Alaska. Whitledge, T. E., and V. A. Luchin (1999), Summary of chemical distributions and dynamics in the Bering Sea, in Dynamics of the Bering Sea, edited by T. R. Loughlin and K. Ohtani, pp. 217 – 250, Univ. of Alaska Sea Grant, Fairbanks. Whitledge, T. E., W. S. Reeburgh, and J. J. Walsh (1986), Seasonal inorganic nitrogen distribution and dynamics in the southeastern Bering Sea, Cont. Shelf Res., 5, 109 – 132.

C. J. Deal, M. Jin, N. Tanaka, and J. Wang, International Arctic Research Center, University of Alaska Fairbanks, Alaska 99775-7340, USA. ([email protected]) M. Ikeda, Graduate School of Environmental Earth Science, Hokkaido University, Sapporo 060-0819, Japan. ([email protected])

20 of 20