AN EVALUATION OF A COMPUTER SIMULATION MODEL OF PLANKTON DYNAMICS IN MONTEREY BAY David Edward Henri ckson NAVAL POSTGRADUATE SCHOOL Monterey, California - THESIS AN EVALUATION OF A COMPUTER SIMULATION MODEL OF PLANKTON DYNAMICS IN MONTEREY BAY by David Edward Henrickson September 1976 Thesis Advisor: E. D. Traganza Approved for public release; distribution unlimited. T176098 SECURITY CLASSIFICATION OF THIS PACE rDltn Dmtm Inrtr.a? REPORT DOCUMENTATION PAGE READ INSTRUCTIONS BEFORE COMPLETING FORM 2. GOVT ACCESSION NO. 3. RECIPIENT'S CATALOG NUMBER 4. TITLE and Sub(llla) An Evaluation of a Computer Simulation Model of Plankton Dynamics in Monterey Bay S. TYPE OF REPORT ft PERIOD COVERED Master's Thesis; September 1976 • PERFORMING ORG. REPORT NUMBER 7. AUTHOn(t) David Edward Henrickson ft. CONTRACT OR GRANT NUMBERS » PERFORMING ORGANIZATION NAME ANO ADDRESS Naval Postgraduate School Monterey, California 9 3940 10. PROGRAM ELEMENT. PROJECT, TASK AREA ft WORK UNIT NUMBERS II. CONTROLLING OFFICE NAME ANO ADDRESS Naval Postgraduate School Monterey, California 93940 12. REPORT DATE September 1976 IS. NUMBER OF PAGES 93 14. MONITORING AGENCY NAME ft ADDRESS*"// dttlaront from Controlling OHIetj Naval Postgraduate School Monterey, California 93940 IS. SECURITY CLASS, (ol thla rdpert) Unclassified IS*. OECL ASSIFI CATION/' DOWN GRADING SCHEDULE 16. DISTRIBUTION STATEMENT (ol thla Roporl) Approved for public release; distribution unlimited. 17. DISTRIBUTION STATEMENT (oj iha omatrmct tncrmd In Block 20, II dlllmrmnl fooai Ropori) '8. SUPPLEMENTARY NOTES U. »GErn^(n Dnm Erti»r»d based on the generation of both advective fog and low stratus cloud cover common during upwelling on the California coast. Analysis of the model's response to sinking and advection of phytoplankton was examined. The importance of seasonal in- creases in predators was introduced as a controlling factor in the seasonal growth of zooplankton. The model is able to predict the seasonal trends of phosphate, phytoplankton, and zooplankton throughout the year. J DD 1473 Form 1 Jan 73 S/N 0102-014-6601 2 SECURITY CLASSIFICATION OF THIS PAGEf***" Dalm Enltrmd) An Evaluation of a Computer Simulation Model of Plankton Dynamics in Monterey Bay by David Edward ^enrickson Lieutenant, United States Coast Guard B.S., United States Coast Guard Academy, 1971 Submitted in partial fulfillment of the requirements for the degree of MASTER OE SCIENCE IN OCEANOGRAPHY from the NAVAL POSTGRADUATE SCHOOL September 1976 H45-H DUDLEY KNOX LIBRARY HAVAL POSTGRADUATE SCHOOL MONTEREY. CALIFORNIA S3S40 ABSTRACT A computer simulation model of the phosphate, phytoplank- ton and zooplankton dynamics in Monterey Bay was examined and modified. The model is driven by four forcing functions ex- pressed as annual cycles of upwelling velocity, incident solar radiation, mixed layer depth, and mixed layer tempera- ture. An alternate upwelling index was developed based on the local wind field. A revised radiation index is employed based on the generation of both advective fog and low stratus cloud cover common during upwelling on the California coast. Analysis of the model's response to sinking and advection of phytoplankton was examined. The importance of seasonal in- creases in predators was introduced as a controlling factor in the seasonal growth of zooplankton. The model is able to predict the seasonal trends of phosphate, phytoplankton, and zooplankton throughout the year. TABLE OF- CONTENTS I. INTRODUCTION ------------------ 10 A. PURPOSE ------------------ 10 B. MODELING PHILOSOPHY IN THE ECO- LOGICAL CONTEXT -------------- 10 C. DESCRIPTION OF THE MODEL ---------- 13 II. BACKGROUND THEORY --------------- 16 A. UPWELLING ----------------- 16 1. Ekman Dynamics ------------- 16 2. Wind Stress on the Sea Surface ----- 18 3. Spatial Extent ------------- 19 B. CURRENTLY AVAILABLE UPWELLING INDEX - - - - 2 2 C. EFFECTS OF UPWELLING ON INCIDENT RADIATION ----------------- 22 D. SINKING OF ALGAL CELLS ----------- 24 E. ADVECTION ----------------- 26 F. PREDATION PRESSURE ------------- 29 III. METHODS -------------------- 33 A. UPWELLING PROGRAM ------------- 33 B. RADIATION INDEX -------------- 35 C. ALGAL SINKING TERM ------------- 38 D. PHYTOPLANKTON ADVECTION TERM -------- 39 E. PREDATION LOSS --------------- 41 IV. RESULTS -------------------- 46 A. UPWELLING CALCULATIONS ----------- 46 B. SENSITIVITY ANALYSIS ------------56 C. COMPARISON OF SIMULATION RESULTS WITH OBSERVED DATA -------------70 V. DISCUSSION -------------------75 VI. SUGGESTIONS FOR FURTHER RESEARCH --------77 APPENDIX A APPENDIX B APPENDIX C Monterey Bay Upwelling Program ----- 79 Monterey Bay Ecosystem Model Program - - 83 Forcing Functions Used in Simulation - - 86 LIST OF REFERENCES ------------------90 INITIAL DISTRIBUTION LIST --------------92 LIST OF FIGURES Figure 1. Pictorial representation of Monterey Upwelling Ecosystem ---------------- 1M- 2. Upwelling index for 19 7 4 at 36 degrees North, 12 2 degrees West, computed by Bakun - - - - 23 3 3. Monthly means of upwelling index (m /sec) and percent of possible sunshine (after Tont , 1975) ----------------25 4. Predation pressure as a function of prey density (zooplankton biomass) (after Patten, 1971) ---------------31 5. Percent of possible sunlight as a function of upwelling index .----------------36 6. Seasonal variation of incident solar radia- tion, clear sky conditions (Q ) and reduced by fog and clouds (Q.) --------------37 7. Predation pressure on zooplankton by transient predators ----------------42 8. Predator population ----------------44 9. Frequency distribution of wind field direction at Monterey Peninsula Airport , 1974. Wind direction was corrected for topographic effects ----------------47 10. Upwelling index for Monterey from cor- rected airport wind observations, 1974 ------ 48 11. Seasonal temperature variation off Monterey Bay. Values are the mean of four to nine stations (Traganza et al., 1976) ----- 49 12. Seasonal phosphate variation off Monterey Bay. Values are the mean of four to nine stations (Traganza et al . , 1976) ---------51 13. Seasonal salinity variation off Monterey Bay. Values are the mean of four to nine stations (Traganza et al . , 1976) ---------52 14. Chart of the Monterey area showing station location -----------------54 15. Possible "fountainhead" effect over canyon head --------------------55 16. Simulated seasonal phosphate concentration with various phytoplankton sinking rates ; no predation -----_-----________57 17. Simulated seasonal phytoplankton concentration with various phytoplankton sinking rates; no predation -------------------58 18. Simulated seasonal zooplankton concentration with various phytoplankton sinking rates; no predation ------------------- 59 19. Simulated seasonal phosphate concentration with various phytoplankton sinking rates; with predation ------------------61 20. Simulated seasonal phytoplankton concentration with various phytoplankton sinking rates; with predation ------------------62 21. Simulated seasonal zooplankton concentration with various phytoplankton sinking rates; with predation ------------------63 22. Simulated seasonal phosphate concentration with various phytoplankton advection rates (with predation and algal sinking rate = 1.0 m/day) __--_____-_--------- 65 23. Simulated seasonal phytoplankton concentration with various phytoplankton advection rates (with predation and algal sinking rate = 1.0 m/day) -------------------- 66 2M-. Simulated seasonal zooplankton concentration with various phytoplankton advection rates (with predation and algal sinking rate = 1.0 m/day) --------------------67 25. Simulated zooplankton response to various pre- dation conditions (algal sinking rate = 1.0 m/day and advection coefficient = 0.1 m"-1-) - - - - 69 25. Comparison of simulated seasonal phosphate with observed data ----------------71 27. Simulated seasonal phytoplankton concentration - - 72 28. Comparison of simulated seasonal zooplankton concentration with observed data ---------73 ACKNOWLEDGEMENTS The author wishes to thank Dr. Eugene D. Traganza for allowing the author to join in his research project: Inves- tigation of 3iochemical Relationships for Determining Concen- tration of Zooplankton Biomass and its Correlation with Chemical and Acoustical Properties of the Ocean. This pro- ject was supported by the Office of Naval Research, Code 4-82, Arlington, Virginia. Sincere appreciation is extended to Dr. Glenn Jung and Dr. Alden Chace for their helpful comments and pertinent re- commendations in reviewing this research. The author wishes to thank Dr. Robert Renard and the meteorology staff of the Naval Postgraduate School for assist- ance with obtaining and reducing wind data. Finally, the author's sincerest gratitude is felt for his wife, Nancy, whose personal sacrifices during the past nine months have made completion of this project possible. I. INTRODUCTION A. PURPOSE In the interests of assessing the biomass of oceanic areas, various mathematical models have been designed to predict the dynamic response of ecosystems to change in the physical and chemical environment. Such models necessarily rely on accurate characterization of the forcing functions and accurate representation of the system with equations which are consistent with the specific region under study. Recent attention paid to modeling systems in upwelling regions, Coastal Upwelling Ecosystems Analysis (CUEA),has brought about developments in both these areas. Adequate data of a physical and chemical nature is becoming available and significant progress in the refinement of governing equa- tions is evident (Walsh, 1973). The research presented here was aimed at creating a simu- lation model to describe the seasonal plankton dynamics in Monterey Bay, California, as known from the best available data. An existing simulation model (Pearson, 1975) was evaluated in an effort to make refinements and to judge the applicability of the time simulation technique to the Monterey region, characterized by seasonal upwelling. 3. MODELING PHILOSOPHY IN THE ECOLOGICAL CONTEXT The merit of mathematically modeling systems of a purely physical nature in which parameters, initial conditions, 10 boundary conditions and the variable relationships are readi- ly specified is unique. The value of current models of a biological nature requires some comment. According to Odum (1971) there are a number of reasons behind constructing any mathematical model. Prediction of future states of the variables involved is an end in itself. Use of the prediction accuracy to focus data acquisition efforts or direct attention to deficiencies in concepts is likewise valid. If the model is real in the sense that it attempts to predict changes in the state variables based on a framework drawn from research in the real world, then failures may be useful in delineating weaknesses in the governing equations. Construction of a "real" model is at best difficult. A perfect model of an ecosystem assumes that "all states and rates of change of the variables are known at all times" (Walsh, 1973). Given the vastness of complex interactions in even the simplest of biological systems, it may be stated that "no perfect representation of the real world exists except the real world itself" (Walsh, 1973). It follows then that certain assumptions, logical state- ments and approximations of real dynamics are necessary to allow the model to operate. The resulting model may bear little correlation to the actual physical, chemical and biological situation at hand. Again, the argument for con- structing the model must be considered. For purposes of prediction, a model which makes a logical statement relating 11 two state variables in the attempt to specify a relationship which may be unknown in part or in whole has merit. The fact that the logic of the model is a crude simplification is incidental if the prediction capability is proven. All models can be characterized on the basis of realism, precision, and generality (Odum, 1971). It has been suggested that these three qualities are mutually exclusive in ecologi- cal models (Patten, 1971). A fourth category, simplicity, might be considered. Exacting detail, while lending to realism and precision often forces the model to be specific. Only when the model becomes an isomorph, i.e., a one-to-one correspondence to the real world, will generality be restored (Walsh, 1973). The problems in creating a realistic model have been dis- cussed briefly. The value of this type of model lies in re- search guidance. Where a specific question is asked of an ecological model, precision in prediction capability may be enhanced by sacrificing the other qualities. Such an approach might be taken in a fisheries model where a precise output for a limited region is desired (Odum, 1971). Simplic- ity and generality at the expense of both precision and realism may enhance understanding of broad ecological con- cepts . The study presented herein was conducted with the idea of creating a model having good predictive capabilities for conditions occuring over a long span of time in a limited oceanic region. 12 C. DESCRIPTION OF THE MODEL The original ecosystem model was developed by Pearson (1975) as part of continuing research into the correlation of zooplankton biomass with the chemical and acoustic pro- perties of the ocean (Traganza, 1976). It is the intent of this research to identify areas of weakness in the model and make the appropriate changes to the equations and input parameters. These changes include a refinement of the up- welling index and radiation index as well as inclusion of advection, sinking and predation terms in the governing equations . The modeling technique used is a time dependent simula- tion solved digitally on an IBM 360 computer using the IBM Continuous Simulation Modeling Program (CSMP-360). The dynamic equations are written in non-linear differential form and driven by four exogenous forcing functions expressed as annual cycles of upwelling velocity, mixed layer depth, mixed layer temperature, and incident solar radiation. The output consists of annual cycles of phosphate concentration Cmicrogram-atoms of phosphorous per liter, ug-atP/£), phyto- 2 plankton biomass (grams of carbon per square meter, gC/m ) 2 and standing stock of herbivorous zooplankton (g C/m ) in the mixed layer. The Pearson version (Figure 1) of the model is defined by the following basic equations: 1) IP- = NUT + REGEN - UPTAK (1) dt 13 LIGHT Figure 1. Pictoral representation of Monterey Upwelling Ecosystem. 14 whe^e * dXl -^r— = time rate of change of phosphate concentra- tion (ug-at P/l) in the mixed layer NUT = input of phosphate from upwelling and mixing from below the mixed layer (yg-at P/l) REGEN = phosphate recycled as biological excretory products within the mixed layer (yg-at P/l) UPTAK = phosphate depleted by phytoplankton utiliza- tion (yg-an P/l) 2) =¥- = PROD - RESP2 - GRAZ (2) dt dX2 where: -7- — = time rate of change of phytoplankton biomass dt (gC/m2) 2 PROD = photosynthetic production (gC/m ) 2 RESP2 = respiration losses (gC/m ) 2 GRAZ =' grazing losses (gC/m ) 3) 5S-1 = GRAZ - RESP3 - VOID - LOSS (3) dt GRAZ = input9due to ingestion of phytoplankton CgC/nT) 2 RESP3 = respiration losses (gC/m ) 2 VOID = excretory losses (gC/m ) 2 LOSS = stock depletion by mortality (gC/m ) (Pearson, 1975) 15 II. BACKGROUND THEORY A. UPWELLING THEORY 1. Ekman Dynamics Seasonal upwelling of nutrient-rich deep water into the productive mixed layer is an important aspect of the •Monterey Bay ecosystem. The principal source of phosphate in the mixed layer is vertical advection which is produced in response to wind-induced Ekman circulation. Theoretical calculations of the offshore component of horizontal current are based on Ekman pure drift theory in an infinitely deep, homogeneous ocean (Neumann and Pierson, 1966). The compo- nents of the horizontal current take the form: U = V. e-(7T/D)z C0S[45-U/D)Z] and (4a) o V = V e"(TT/D)Z SIN[45-(tt/D)Z] (4b) o D = (A/pu) SIN <|>)'s and (5a) V = T /Aa (5b) o y where U = velocity component in the x direction V = velocity component in the y direction (parallel to wind) V = surface velocity component (45 to right of wind j Northern Hemisphere) Z = depth D = Ekman depth of frictional resistance A = coefficient of eddy viscosity 16 t = surface wind stress in y direction (dynes) a = fp/A

The coefficient of eddy viscosity, A, was assumed constant since detailed information on the small scale mo- tions of the surface layer was not available. This assump- tion is generally made in mass transport studies (Neumann, 1968). The mass transport due to wind driven current is found by integrating the equations of velocity (U and V) over the depth of the water column, as shown: S = p UdZ = t /f (6a) x y and S = p VdZ = 0 (6b) y where : S = mass transport in the x direction S = mass transport in the y direction From the equations for mass transport, it is observed that the net transport is 90 degrees to the right of the wind in the northern hemisphere and that it is directly proportional to the wind stress acting on the sea surface, parallel to the coast. 17 2 . Wind Stress on the Sea Surface The coupling of wind energy to the water at the air- sea interface is defined by the conventionally-accepted form of the wind stress equation : T = PfCzW^ (Gordon, 1972) (7) 2 where: T = wind stress parallel to the coast (dynes/cm ) p1 = density of air C = drag coefficient at height z (non-dimensional) z W = wind speed at height z (cm/sec) As in other momentum exchange studies , the dynamic coefficients appear to present the singular most difficult problem. Wilson (1960) states that the drag coefficient for air flow over a water surface is dependent on the wind speed (W) , the height of measurement (z ) , a surface roughness para- meter (z ), and atmospheric stability terms. There may be additional dependence on oceanic parameters of depth, fetch distance and wave height. An average value of the drag co- efficient for high wind speeds (i.e. greater than 10 m/sec) -3 of 2.37 X 10 is arrived at by examination of research com- pleted through 1959 (Wilson, 1960). Work by Deacon (1962) established an empirical rela- tionship for winds up to 13 meters/second which gives a linear dependence of C (at z = 10m) to the wind speed as follows : C1Q = (1.10 + 0.04 W ) X 10"3 for W in m/sec (8) (Roll, 1965) 18 Deacon's research was carried out using data from ship obser- vations and coastal regions. 3. Spatial Extent The offshore movement of water in a region bounded by a coastline causes replacement water to be upwelled from below the layer affected by the surface wind stress. Accord- ingly, it is necessary to specify the spatial dimensions of the region in order to apply principles of continuity in de- termining the vertical current speed. The significant para- meters involved in upwelling are: (1) the offshore horizontal dimension of the upwelling region, which is most directly related to the width of the wind field; (2) the offshore component of the Ekman current (U) ; and (3) the depth of the Ekman layer (D) . Estimates of the maximum depth of the source of up- welled water in a coastal region are approximately 100 to 200 meters (Neuman and Pierson, 1966) based on the slope of the isotherms. The mass transport offshore due to wind drift current has been previously detailed as occurring in a surface layer of depth, D, corresponding to the Ekman depth of frictional resistance; regardless of the source depth of the upwelled water, it transgresses the "boundary" at depth, D, before being carried offshore. The vertical dimension of the region is then readily specified. It is assumed that the depth of frictional resistance (D) and the depth of the mixed layer (Z) are approximately coincident. The differ- ences between (D) and (Z) are important when phytoplankton 19 are advected out of the Ekman layer but exist throughout the mixed layer. The horizontal extent of the upwelling region is considerably more nebulous than the depth of the mixed layer. Classical estimates of this dimension limit the zone to probably no more than 10 0 kilometers in offshore width. It may be expected that the width of the region is dependent on several factors; among them, the width of the wind field, variations and non-linearities in the energy exchange pro- cesses due to surface roughness characteristics or stability changes in both the atmosphere and the water column, as well as spatial and temporal oscillations of the significant wind vector. In addition, local coastal and bottom topography may figure extensively in the problem as will be discussed later. Sverdrup (1938) observed a relatively well defined offshore boundary to a coastal upwelling zone, coincident witha downwind current which was marked by an intense verti- cal gradient of velocity (Sverdrup et al., 19M-2). He further observed an offshore migration of the current band at a rate somewhat less than the speed of the offshore surface current within the upwelling region. The latter fact implied the possibility of cellular circulation patterns in the near sur- face upwelling zone. Hidaka (195M-) proposed a steady state theory of up- welling in which he defined a horizontal frictional distance, D, , analogous to the depth of frictional resistance 20 appearing in Ekman's work. Applying reasonable values to Hidaka's expression confines the total cellular circulation region to a width of about 15 kilometers. The width of the region is defined by: DR = tt(2Ah/2o) SIN )2 where: AR = 107(gcm_1 sec'1) (9) (Smith, 1968) Downward vertical currents occurring in the seaward half of the cell limit the upwelling to a width of 7.5 kilometers. The theoretical width is not consistent with observations of upwelling at greater distances from the coast (Barnes, 1969). Yoshida (1967) developed a quasi-steady state model applicable to eastern boundary current regions (Smith, 19 6 8^. An expression was derived for the horizontal dimension of the coastal upwelling region which is given by: L = (gH^p/p) (Hy/Dy)^ for latitudes greater than x x 22 degrees (10) where: L = horizontal width of the coastal boundary region (m) _2 g = gravitational acceleration (m-sec ) H = thickness of the upper layer (m) D = H + thickness of the lower layer (m) f = Coriolis parameter (sec" ) _3 p = density of water (g-cm ) Ap = density difference between deep and surface layers (g-cm-3) y = internal friction coefficient (sec ) Y = 2 X 10"7 sec"1 (Smith, 1968) 21 3. CURRENTLY AVAILABLE UPWELLING INDEX Bakun (197 6) has calculated coastal upwelling indices for the coast of North America using Fleet Numerical Weather Central monthly mean pressure fields to compute the geos- trophic wind. Analysis is done on a 63 x 63 point grid of approximately 200 nautical mile mesh length (Bakun, 1973). Figure 2 shows Bakun' s results for 1974 at 36 degrees North, 122 degrees West, approximatley 5 0 miles south of Monterey Bay. In reviewing several aspects of hydrographic surveys taken in Monterey Bay, it was suspected that this region is characterized by local upwelling patterns not appearing in the index computed by Bakun. Further discussion follows in the methods section. C. EFFECTS OF UPWELLING ON INCIDENT RADIATION Cold water upwelled from depth during the spring and summer months brings about both advection fog and stratus cloud formation. The fog occurs as warm surface air is cooled by the cold seawater to its saturation point causing the condensation of the contained water vapor. Stratus clouds occur below the base of a quasi-permanent atmospheric temperature inversion (Tont, 197 5). The net effect of both fog and stratus layers is to re- duce the amount of solar radiation reaching the sea surface during the upwelling period. Tont (1975) obtained mean values of solar radiation at the surface over the period 1950-1973 at San Diego. A correlation study of the annual upwelling index and the amount of sunlight at the surface 22 -M U o CO CD CD in bO CD U3 CD C CQ >^ rQ 4-> n3 X) CD =ir -H 3 Oh e o o en U O 4h CD C •H ■H CO CD IS 03 CD bO CD C & bO H CD rH TD 0) D H CM 0) h bO ■H (peieDS) X3QN I 9Nlll3Mdn 23 showed a maximum of about M-0 per cent reduction in available sunlight during the peak upwelling months. The result of Tont ' s research is shown in Figure 3. D. SINKING OF ALGAL CELLS Review of sedimentation rates of algal cells indicate justification for including a sinking term in the phytoplank- ton equation. The Pearson (1975) model assumes uniform algal concentation within the mixed layer. It is reasonable that changes in the concentration of phytoplankton caused by ver- tical movement out of the layer will occur. Parsons and Takahashi (1973) state that the theoretical sinking rate of phytoplankton can be determined from: v = ifli (£l=£) (ID where : r = radius of the cell g = gravitational acceleration p' = density of the organism p = density of the medium (sea water) r\ - viscosity of the medium = form resistance coefficient (accounts for non- r spherical shapes) Measured rates of sinking for live algal cells vary be- tween zero and 3 0 meters per day according to Parsons and Takahashi (1973). These rates represent a broad range of soecies and sizes. 24 3N1HSNPIS 3"l9ISS0d lN30d3d o CD & CO CD -H CO 4h £ «d v^ s_^ X CD CD C T3 •H C rC ■H bO 3 C w •H H CD H H CD £1 3 •H ft W 3 en 0 Uh ft o 4-1 01 0 c • (d ■M /— \ CD C LO 6 CD r- 0 CD >, f- H CD ,C ft *N -M -P c -n c 0 C 0 s rd H • 00 CD h 3 b£ •H Fm O o o o i-oas-puj) X3QN o o 9NI"113Mdn 25 E. ADVECTION It is known that significant current velocities occur within the mixed layer from Ekman dynamics. It might be ex- pected then that the concentration of phosphate, phytoplank- ton and zcoplankton will be affected by advection. O'Brien and Wrobleski (1972) have shown by scale analysis that advec- tion and biological productivity of the Peru upwelling eco- system are equally important. The advective terms of the material derivative of a property (C) are written: Advective change =u|§+v|y+WI§ (12) where: U,V,W = velocity components in x, y, and z directions respectively 3C 3C 3C ttt, iry, -yY - gradients of the concentration of property (C) in the x, y, and z directions A review of California Cooperative Fisheries Investiga- tion (CALCOFI) data for 1974 (CALCOFI, 1974) shows that the phosphate concentration varies significantly in the vertical direction but the horizontal variation is not as significant. It is assumed for the purposes of this model that the ad- vection term for phosphate reduces to the vertical component: 8X1 Advective change = W js— (13) 3X1 3X1 n since: ^- = ^- - 0 Pearson (1975) has included this term in the computer model as: UPWELL = W DX1 ~ X1 (14) Jl J\J_i 26 where : W = vertical speed (m/day) EKL = depth of the Ekman layer (m) , EKL = D DX1 = concentration of phosphate below the mixed layer (yg-at P/l) XI = concentration of phosphate in the mixed layer (constant) (yg-at P/l) DX1 - XI gY? = vertical gradient of phosphate (yg-at P/m 1) The advection of zooplankton in the revised model is discounted because it is assumed that they have developed mechanisms which permit them to remain in a particular oceanic region by swimming or riding cellular circulation currents . Phytoplankton, however, are affected by advection. CALCOFI data (CALCOFI, 1974) shows that the significant gradients of phytoplankton occur in the vertical and off- shore directions; therefore, the advection expression for phytoplankton is reduced to: Advective change = U rrrr— + W -^y— (15) The horizontal velocity is related to the vertical veloc- dV dY ity by continuity (when it is assumed that -jy- = 0, parallel to the coast) such that: U = - WL/D (16) 27 where: U = horizontal velocity component offshore (m/day) W = upward vertical velocity component (m/day) L = width of the region (m) D = Ekman depth (m) The phytoplankton advection term can be written: Advective change = W [- | ||^ + ||i] (17) The gradients of phytoplankton can be approximated by: 9X2 LIM AX2 , 3X2 _ LIM AX2 MQ ... 9X~ ' AX+0 AX~ ana 3Z~ " AZ+0 AZ~ Qiaafet; If it is assumed that the concentrations of phytoplankton below the mixed layer and outside the modeled area are much less than the concentration in the modeled region, then the gradient can be approximated by: = horizontal gradient; (19) AX X2-0 (20) ■j-j— - vertical gradient; where: X2 = phytoplankton concentration in the mixed layer AX = horizontal distance at which the concen- tration of phytoplankton becomes negli- gible with respect to the modeled area AZ = vertical distance at which phytoplankton concentration becomes negligibly small The phytoplankton advection term can then be written as : Advective change to , , Phytoplankton in the = W(X2)[- ~% + aY] (21) modeled area 28 The sign convention is established as positive (+) into the model area and negative (-) out. F. PREDATION PRESSURE Under constant environmental conditions , the dynamics of zooplankton growth are a function of the food supply and predation pressure. Natural processes of population control insure that a given species will not be eliminated by low food supply and high predation pressure or that unstable growth will not result from the inverse situation. Attempts to translate the stabilizing factors into mathematical terms often result in simplified general relationships that do not exhibit the flexibilities inherent in the real system. Not- withstanding the limitations, a predation term is included in the simulation model. The long term survivability of both predator and prey is keyed directly to maintenance of a balance in energy expendi- tures and gains. To date, laboratory experiments to dupli- cate this balance have not been entirely successful. Prey stocks in small laboratory environments have been artificial- ly supported by providing refuge where predator access to the prey was denied and by introducing additional prey to the system to replenish the supply (Patten, 1971). Neverthe- less, population control by predation is prey-density depen- dent as a first approximation. For a system of one prey species and a single predator species, predation pressure may be represented by: P = k1 D for 0d1 (22) 29 where: P = predation pressure or food intake per unit predator D = prey density d = satiation threshold of prey density k, and k« = constants Above prey density values of d, , predator satiation will oc- cur and pressure will level off at some constant high value. The function is shown on an arbitrary scale in Figure 4 according to Patten (1971). This simple relationship is not entirely satisfactory for the predator since low food supplies would result in starvation. Under such circum- stances, and given a second prey species of perhaps less palatability , the predator would switch food intake to the more abundant species of prey. There are several obvious implications to the switching phenomenon: 1. predation on the most abundant species assures stable growth of that species , 2. switching from a declining species serves to prevent elimination of that species (Patten, 1971), and 3. elimination of the predator by overgrazing one avail- able food source is avoided. The model under consideration in this thesis is a single species approximation in which predation pressure is incor- porated with other zooplankton mortality factors, i.e. morphological death, in the "LOSS" term. This term is given by Pearson (1975) as: LOSS = L (X3) (23) 30 PREY DENSITY Figure U- . Predation pressure as a function of prey density (zooplankton biomass) (after Patten, 1971). 31 where: LOSS = fraction of the zooplankton population re- moved from the system in a day's time (g C/m2 day) L = constant rate of loss or percent loss per day (day~l) o X3 = zooplankton "standing crop" (g C/m ) Harriston concluded in 1960 that "herbivores are preda- tor limited" (Patten, 1971). Conclusions by Patten (1971) show that herbivores are both food and predator limited. Therefore, in the modified model, the LOSS term is reserved for variable predator limi- tation of herbivorous zooplankton and natural death of zooplankton herbivores is assumed to be negligible. 32 III. METHODS The existing model appears weak in several areas. Orig- inally, the forcing functions of upwelling and solar radia- tion which are critical to the simulation were obtained from average cycles for the California coast (upwelling) and the latitude band between 3 0 and M-0 degrees north. Phytoplankton sinking and advection terms were lacking and the predation pressure term in the zooplankton equation needed improvement. A. UPWELLING PROGRAM A computer program was written to calculate the annual upwelling index for Monterey Bay. Wind data was obtained from the Monterey Peninsula Airport observations during 1974. Hourly observations were examined to determine the signifi- cant mean daily wind vector for each day and corrections were applied to speed and direction when indicated by topo- graphic obstructions. Fortunately, no corrections were needed to northerly winds which are the driving mechanism for upwelling. An x-y coordinate system was established such that the 2 surface wind stress (x = p'C W ) could be directly calcu- y lated from the wind component parallel to the coast. The Ekman mass transport (S ) , surface current velocity (V ) and the average velocity within the Ekman layer (U) were calculated by the following relations : V = x /Aa (5b) o y 33 where: V = surface current vector o and ■»/ U = ± / UdZ (24) where: U = average velocity in x direction in the layer and: S = x /f (6a) x y v ' where S = mass transport in x direction (offshore) The vertical current representing upwelling is determined per unit y (along coast) from the continuity equation, (7 - v) = 0, such that 0(D) = W(L) (25) where: U .= average horizontal velocity in the Ekman layer (m/day) W = upwelling speed (m/day) D = Ekman depth (m) L = horizontal width of the upwelling region (m) Zero current is assumed parallel to the coast. The horizontal width of the upwelling region, (L) was determined from Yoshida's (1967) equation. Applying typical values to this equation yields a horizontal width of 3 X 10 m. Principles of continuity were used to arrive at the mean daily upwelling current values and these figures were aver- aged over a 3 0-day period. 34 B. RADIATION INDEX A revised incident solar radiation index which takes in- to account the reduction of sunlight by local fog and stratus cloud cover was developed. Correlation of the measurements of surface irradiance and the strength of upwelling (as in- dicated by the upwelling index) by Tont (1975) is repre- sented by Figure 5. This figure is derived from the results of Tont ' s study as shown in Figure 3 by plotting the percent of possible sunshine (Y-axis) as a function of the upwelling index (X-axis) which has been scaled to a maximum value of one. The resulting curves (A) and (B) show the conditions occurring after, and prior to, the upwelling maximum, respec- tively. The difference exhibited between the curves may be due to changes in the character of the air mass brought about by seasonal migrations of the quasi-permanent thermal low over Nevada and the North Pacific sub-tropical high located west of the coast. Figure 5 was used to calculate a revised incident radia- tion index by entering with the newly-developed upwelling index values and by applying the corresponding percent to the theoretical mean monthly radiation for clear sky condi- tions (possible sunlight). Figure 6 shows the possible sun- light throughout the year (from Tont, 197 5) as Q , and the theoretical radiation at the surface (Q.) after the effects l of fog and clouds are considered. 35 0.2 0.4 0.6 UPWELLING INDEX 1.0 Figure 5. Percent of possible sunlight as a function of upwelling index. 36 o >. X CO in rd CD • !~i *•■> 0 •H cy «%_• c 0 CO •H X3 +J 3 iK o •H rH T3 0 03 fc T G fc rd rd ^ c ,q CD T3 T3 •H a) a 0 c ^ •H T co e H 0 id ■H c ■P o •H CO T3 rd c CD 0 GO o CO CD W3 37 The incident solar radiation is an essential part of the ecosystem simulation model in that it is used in the equations governing the production of organic carbon by the photosynthetic activity of phytoplankton. These equations are discussed in detail by Pearson (1975). C. ALGAL SINKING TERM The phytoplankton sinking term which was added to the Pearson (1975) model is written in the form used by Riley (1965) in the units of the quantity of phytoplankton trans- 2 ferred per day (g C/m -day). The amount of phytoplankton that sinks out of the mixed layer in a day is given by: SINK = (V/Z)X2 (26) where: SINK = flux of phytoplankton out of mixed layer (gC/m2-day) V = sinking rate (m/day) Z = mixed layer depth (m) 2 X2 = phytoplankton biomass in mixed layer (gC/m ) V/Z = fraction of phytoplankton which sinks out of mixed layer per day Vertical circulation with the water column becomes significant as the upwelling current speed and sinking rates of phytoplankton approach the same order of magnitude. Since the phytoplankton are non-mobile and are generally of the same density as the water or have developed shapes which increase their sinking resistance, they will be carried along with vertical currents in the water column. The vertical circulation during upwelling opposes sinking but may 38 accelerate downward transfer during brief periods of surface convergence and associated downwelling. The complete phytoplankton sinking term is: SINK = ((V - W)/Z)X2 (27) where: W = upwelling speed (m/day) An average value of one meter per day was used "for the sink- ing rate (V) in the computer simulation model. This value was determined from estimates presented by Lehman et al . (1975) and Bannister (1974). The SINK term acts to decrease the phytoplankton concentration in the modeled region and is subtracted in the phytoplankton equation (2) as shown: ^-rP- = PROD - RESP2 - GRAZ - SINK (28) dt D. PHYTOPLANKTON ADVECTION TERM An equation describing the horizontal advection of phyto- plankton is included in the revised model. This term is in the units of flux and is written as : ADVEC2 = W (X2)K (29) where: ADVEC2 = the change of concentration of phytoplank- ton over time (in the mixed layer) due to advection (g C/m^-day) W = upwelling speed (m/day) X2 = phytoplankton concentration in the mixed layer (g C/mO K = advection coefficient (m ) 39 The simplified advection term used in the model is de- rived from the equation (21) shown in the background theory- section, which is: Advective change = W (X2) [- ^ + -^] (21) The Ekman current offshore which results in an upwelling current (W) is effective over the depth of frictional resist- ance (D) , but the phytoplankton in the model (which are uniformly distributed in the mixed layer) will be advected at the average velocity over the depth of the mixed layer as determined by the fraction of the mixed layer that is coinci- dent with the Ekman layer. It is known that the mixed layer varies seasonally from about 10 to 100 meters in depth and this causes an annual (seasonal) variation in the average horizontal current in the phytoplankton advection term. The average current of the mixed layer can be determined from: z U = — / U dZ where: U? = a function of e (30) A shallow mixed layer experiences higher average veloci- ties than a deep layer as seen from the integral expression (eq. 30). A coefficient to describe the seasonal average of the vertically averaged horizontal current in the mixed layer is included in the advection equation as shown: Advective change = W (X2) [- j^ Km + ^] (31) 40 The coefficient, K , represents the seasonal average of the fraction of the mixed layer which is coincident with the Ekman layer and is affected by the Ekman current. The constants in the advection expression (within the brackets of eq. 31.) are lumped into the advection coefficient, K, for convenience in equation (29). The range of K varies -1 -2 -1 from about 10 to 10 m depending on the estimated values of the gradient distance, (AX) and (AZ) and the velocity co- efficient, K . A value of 0.1 m was used in this model ' m for K. E. PREDATION LOSS The predation pressure function hypothesized for Monterey Bay is based on the assumption that the system supports a resident population of herbivorous zooplankton predators throughout the year. The pressure exerted on the herbivore prey by this maintenance level (N) of predators follows the curve in Figure 7 below the prey density threshold, d. The pressure is sufficiently low to allow growth of the zooplank- ton stock, but high enough to stabilize growth. A second pressure is imposed on the zooplankton (see Figure 7) after their density reaches the critical level, d. The simulation approximates conditions which may be brought about as tran- sient predators move into the area presumably in response to increased food availability. In summary, predation pressure may be represented by a two part function with a pressure- density curve, a, due to resident predator species and curve b due to the transient predator population. 41 Ld q: Z) uo oo LlI en CL O < Q LJ C£ DL PREY DENSITY Figure 7. Predation pressure on zooplankton by transient predators. 42 Predator biomass may be expressed as shown in Figure 8, as a simple step increase or some other function of prey density. Consideration of the step type function is impor- tant because it allows the predator population (and conse- quent pressure) to maximize during those seasons when the zooplankton "standing crop" will support additional numbers of predators; the high predator pressure will deplete the food stock to a lower level than would be reached by the resident predator pressure alone. Measurements of zooplank- ton biomass for 1974 (Traganza, 1975) suggest a rapid decline in "standing crop" following an early summer maximum. One might suspect that a transient predator population is forced out of the region once the food sources are depleted. This is done in the model by decreasing the predator population when the zooplankton biomass declines to a specified level. The predator levels were related directly to zooplankton "standing crop" by triggering an increase in predator popula- tion at a threshold of zooplankton biomass during periods when this prey population was on the upswing. The predators were similarly reduced at a second threshold during the de- clining phase of zooplankton growth. A set of four condi- tional statements is used in the simulation to describe the predator-prey relationship. These statements specify the following : IP ^1 > o AND X3 < 1.0 THEN: FISH = 1.0 dt __ dX3_ > 0 AND X3 > 1.0 THEN: FISH =5.0 dt 43 Z) CL O Q_ O I— qN LU cn CL 7 / IT / _ / // / // PREY DENSI T Y Figure 3. Predator population. 44 IF =$£ < 0 AND X3 > 0.2 THEN: FISH =5.0 at — — Ayr* IF ^=~- < 0 AND X3 < 0.2 THEN: FISH = 1.0 at — The modified predation pressure term in the zooplankton equation of The model is: LOSS = (L)(FISH)(X3) (32) where LOSS = the amount of herbivorous zooplankton biomass lost per day as a result of predation (g C/m2 day) L = loss rate or percent loss per FISH per day (day-1) FISH = number of predators (non-dimensional) 2 X3 = zooplankton "standing crop" (g C/m ) 45 IV. RESULTS A. UPWELLING CALCULATIONS The mesoscale wind data used to calculate u'pwelling for this simulation is shown by Figure 9. The plot depicts the frequency of occurrence (f = N/SN) for 4-5-degree segments of the compass. The prominent wind direction is as usual from the northwest .during the April to September upwelling period and shifts to the south in the first and last quarters , but this year there was a significantly high frequency of northwesterly winds during the last months of the year. The computer generated upwelling index (Figure 10) shows a pri- mary upwelling maximum in May, a second peak in about mid September and another in mid November. The possibility of a secondary upwelling peak is suggested by temperature obser- vations by Traganza et al . (1976). Figure 11 shows a rise in the isotherms peaking 30 to 45 days after the indicated wind initiated upwelling maximums in May and September. Al- though the water column does not respond instantaneously to surface wind stress, the delay noted here is excessive (Barham, 1957). The delay may be attributable to the assump- tion that the winds measured at Monterey Airport are actually representative of the winds over the bay, when in fact they may not be. The delay may also be caused by the lack of suf- ficient data to accurately depict the seasonal trends of temperature in the water column. 46 >^ u tu 0 U 4-1 CU TD C 0) 0 ■P O a; +J & nfl in O c a 0 ■H oo P m O s 0) fH c ■H 0 T3 •H P T3 a rH cu 0) Sh •H ■H Hh T3 T3 T3 C c ■H •H 2 3: <4H 0 • e l> 0 en ■H rH • ■P to 3 *\ -p ,0 P o ■H H , rd ,C a rH a C 3 tO cu CO H 3 C b0 ct-h O cu C Ph & CU 0 PlH CU ■p CD 0) Fh b0 •H 47 T3 C •H -H u o PU & •H CD -P a cu p. Ph o o 6 o h CU p. ■H P. d rH . id >, rO ■H PQ CD >^ id CD N Fh c •H c 0) Fh 0 3 4-> ■fj id Fh Fh 3 CD 0 P.I-H £ CD 4h -p O rH a rd rd C CD 0 e w Id CD 0) f-J CO •p o o o o o o o o \n o S83JL3W 04 o o CD Fh bO •H 49 Traganza et al. (1976) observed a sharp rise in phosphate (Figure 12) with the major upwelling in May and a slight in- crease in phosphate and salinity (Figure 13) in September. No salinity data were available for the first half of the year. There was no clear-cut correlation of either phosphate or salinity to the upwelling index during November. Nutrient data from four to nine stations taken during seven cruises in 1974 were averaged (see Figure 12) and show a rapid rise in October which is probably not related to upwelling. From the five year study of Monterey Bay by Bolin (19 64) it is known that the nutrient concentration is characteristically low over the depth of the euphotic zone (0-200 meters) for two or three months at the end of the year. The restoring of the nutrient level in all but the upper 20 to 30 meters to half of its May upwelling value while the surface tempera- ture reached apeak, may have signalled the beginning of the Davidson current period (see Bolin and Abbott, 1962). The Davidson current brings a southerly winter oceanic water mass into the Monterey region. This "Davidson water" is characterized by lower surface temperatures and surface salin- ities (due to high amounts of rain) . The deeper water of the euphotic zone has higher salinities and nutrient concen- trations than exist at Monterey. The mechanism which estab- lishes these characteristics may be winter storm mixing oc- curring south of Monterey or upwelling and mixing from below brought about by divergent (cyclonic) eddies formed in the current stream as it moves northward along the coast (Traganza et al., 1976). 50 CO o- CT> H CO *> • 0) rH 3 rd fd • N fd C fd OJ bO rd >! Ch 0) Eh & p OJ n ■H <4H rd ,C Mh U* n w 0 C ,4 rd Ph a; g rH *"■ rd l ,c p* CO o •H P. rd > ■H •H G O •H +J rd p. p c CD ■ c rd O TD a cu P. P (d CO •* O en a-p rd o c G H rd ^ CD G oo -H CO cd C p> O rd P> H X g c g td •H rH CO & (ID cu p. •H Ph /d P-6n 57 >■ < UJ CO 3 0 •H h rd > ,C -P •H 2 c 0 ■H +J id • & c +J 0 c •H ■M o rd C T3 0 0) a fc CX C 0 0 +j c ^ c • r\ (13 CO H (D ft-H O rd +-1 & >. ,C bO O, C •H rH ^ rd c C •H 0 CO CO rd c a; 0 CO -M ^ -a C a) aj ■p rH rd Ph rH 0 3 ■H e >, •H & CO Oh bfl •H ^w/O 6 >- < Q LU CO o •H & rd > 4-> C O •H • ■M C 03 O fc-H +-> -H C rd CU T3 O Q) C fc o o c o -F c o G C o CO It) CD CO . cu rd -H rH rd P. U O O bO N C •H rH ,* C •H CO C O CO +J T3 C cu rd •P rH rd a, H O P -M 6 >^ •H ,C CO o< oo H CU u bO •H P-. zw/0 6 59 The sinking rate was varied in a second test but preda- tor pressure was exerted on the zooplankton by employing the step function predator-prey relationship detailed earli- er. Additional predators were introduced in the simulation o at a zooplankton biomass level of 1.0 g C/m and were removed when zooplankton biomass fell to 0.2 g C/m . Thresh- old values were determined empirically. The zooplankton bio- 2 2 mass maximum is reduced from 1-10 g C/m to 0.1-1.0 g C/m by the addition of predation pressure (Figure 21). The phosphate curves (Figure 19) show a marked response to phyto- plankton growth as evidenced by a lower nutrient minimum following the peak of the summer bloom. A considerable change in the zooplankton response, i.e. lower and earlier, is evidently due to reduction of food sources as phytoplank- ton is allowed to sink out of the mixed layer. The initial effect of increasing the sinking .rate from zero to three m/day on zooplankton is a shift in the occurrence of the peak to approximately 70 days later. Further increase of the sinking rate to six m/day results in a twofold decrease in the carbon biomass of zooplankton with an additional 2 0 day delay of the maximum. It is apparent that the rate of growth of zooplankton is slowed down and the maximum biomass occurs later because of food limitation but that predator pressure is responsible for limiting the magnitude of the zooplankton biomass. The zooplankton peaks also occur earlier when pre- dation is applied and also occur before the phytoplankton peaks (Figure 20). The simulation, therefore requires a positive sinking rate. 60 I o ft CO o •H U > +J •H 2 C o ■H +J id ■H ■ c +j CO rd -H ft CO »•» O co ,C CO ft -H id id c o CO C H 11 M CO CO CO +-> d -p 6 fd •H H CO ft CD CO /d I o - 6 n 61 co J 0 •H P. fd > A o -p •H o 3 CO on concentration ith predation. o ty> >- -P 3 M o < C — no en Q rH 0) Ph-P 0 13 Z. -P P. — LU a c •H X rd C — C -H h- 0 co CO rd C OJ 0 CO -p o T3 C o oj o o ■P (r) concentration wi ith predation. o U) c 2 o >- < 0 ■p •* M co Q 03 -P ^. o 0 bfl UJ •H 3£ rH ^ 03 C — C -H 0 co h- CO 03 C 0) 0 CO -p o ^ o T3 C - < LU 1 /d lD"6n c 0 1 p ^-s ^ ►> c rd id T3 H \ (X S 0 ■P o >i • ■C H & n co 3 CU 0 p> •H rd P. Pi rd > bO G ,c •H p ^ •H c 3 •H CO C 0 rH •H n3 P bO rd H Pi rd P> C T3 cu G 0 rd G 0 G O 0 •H 0) P -p rd ITJ T3 ^ CU Ph Pi co Ph 0 ,c rC Oh •H •H H 3 rd x_*- c 0 CO to CU rd p CU ft) CO Pi X3 G CU 0 p ■H rd P H O 3 CU £ > •H T) C/D rd CM CM (1) P. 3 bO •H En 65 o o en o o < Q LU O O - O c o p M c \ A S Cu a co • G rH 0 •H II su 03 CU > -P +J •H be 3 C ■H C ^ 0 C •H -H - -: nd & iH p aj G 60 - - — - •H - 3 C : •.: -.: L - - CD & ~ - 0) o g a; s > •H T3 CO 03 — & 60 2^/0 P 66 Z^/0 6 c o p X • c ^-n 03 >, rH 03 P.TD 0 \ P e >^ fl o O, • H CO 3 II 0 •H CD & P 03 03 > & rC bfl P C •H •H 3 ^ C G •H o CO ■H P rH oj 03 fc bO P rH G 03 CD O "O c C o 03 0 C c 0 o •H p P X 03 c TD 03 CD H fc ft Oh O o rC N P •H H 2 03 ^_^ C O CO CO CD 03 P 0) 03 CO & TD c CD 0 p •rH 03 P rH a 2 •H T3 oo 03 • CN CD h 3 bfl ■H Ph 67 The predation pressure terms in the zooplankton equation were varied (Figures 25, 26, and 27) while keeping the sink- ing rate at 1.0 m/day and the advection coefficient at K = 0.1m . The conditions of the test were: 1. predation pressure set to zero; 2. predation pressure defined by the linear function of zooplankton biomass: predator pressure = L times zooplankton biomass, where L = 0.01, and the units are g C/m2-day = (day-1)(g C/m2); 3. the linear function in 2, but L = 0.02; 4-. the step function discussed under the methods section of this thesis. In all cases , the resulting zooplankton growth appears food- limited until about day 12 0 when phytoplankton growth begins to show a rate change. Complete removal of predator pressure allows the zooplankton to increase rapidly until food limita- tion (be depletion of phytoplankton) again occurs on day 275. When the predator pressure increases linearly (with L = 0.01), 2 a single zooplankton peak of 3.05 g C/m occurs about day 300. This is not consistent with observed conditions (Traganza, 1976). Doubling the predation pressure rate to L = 0.02 day" severely restricts zooplankton growth but pro- duces a curve with a hint of temporal conformity to the observed zooplankton maximum and minimum biomass levels. When the predator function is triggered at thresholds of zooplankton biomass, the zooplankton exhibit a double peak which is suggested by Traganza' s (1976) data. 68 cu.H ^ '■h CD o o c 4 — (N u ^ O c a; O o il \ — o Q. a; o 00 o o CN O o fT3 c ^ •H | 8° o ■H +-> rd c o •H +J o S +" > CO CD c o c rd id 0) -M rd OT CM CO O 00 o o CN LO O 0) 3^J /o 6 69 C. COMPARISON OF SIMULATION RESULTS WITH OBSERVED DATA The computer simulation was run with the new upwelling and radiation indices and the step function predation pres- sure term detailed in the preceding section. The sinking rate was set at 1.6 m/day and the advection coefficient set to K = 0.1 m"1. A comparison of the simulation results with data observed by Traganza et al . (1976) is shown in Figures 26, 27, and 28. The calculated nutrient .values lag the averaged mixed layer concentrations by 20 to 30 days during the summer months with a moderate amount of error in magnitude. The nutrient maxi- mum occurs on day 144 with a value of 2.08 yg-at P/l in the simulation. A rapid decrease in phosphate levels coincident with the summer phytoplankton bloom displays the impact of nutrient uptake by phytoplankton in the model. The simulated phytoplankton peak occurs 40 days after the nutrient peak during the period of maximum incident radia- tion (see Appendix C). The modeled response of phytoplankton was relatively inflexible with regard to the radiation index, i.e. the phytoplankton peak was found to occur within a few days of the time of maximum incident radiation when the sink- ing rate was varied from about 1.0 m/day to 1.8 m/day. The relationship of the phytoplankton peak to the radiation peak is explained by examination of the phytoplankton growth equa- tions (shown in Appendix B). A slight increase in phytoplank- ton biomass was simulated during December (day 330) which could be due to the increased phosphate concentration at the end of the year. 70 o o 00 o o L0 >- < Q L±J O O O /d }D-6n td +-> T3 •O CD > a> co O £ +J •H CD ■M nj ,C ft CO o ,c a c o CO 03 CU CO CU +J H S •H CO Uh O c o CO •H Ph 03 ft £ o o to CM bO •H 71 o o CO o >- o < LlI o o o •H -H & c cu a c o o c o ■p X c rH a, o >1 x: a a) c o w it( 0) w 0) -H (TJ rH 6 •H CO (1) ■H Pj-t ^lu/0 6 72 o o 00 o o CO >- < Q UJ o o +J •H c o •H ■H id h ■P C 0) o C 0 o c o c id H Oh 0 o N rH n3 C o CO id CO -a a •«■ * 4*- 45- * 4t #■ :* < -* # # # •ft***-*-**-***-* *.«-**^^*# < *--»-*4J-4i--«-4S-4*-4*--»-4!--a-4*-41-##-*.K- 4**4*-4***#4i-4$-* #####*### * * *■* 4* 4*4*4* 4** 4** 4* 4*4* * 4t*# 4* 4* 4*4* 4* *4<- * ■»• 4*4*4* 4* 4*- •* * * 4* 4*-* 4* *- 4fr 4* * 4* 4*4*4* 4t# 4*-* ■»■■»•■«■•»• #*•#-* Z Z *** * # 4** # xi-- ■«■■«■■»•* 4t * ## h>*-*- •»■■«■•«•*• Vr * ** >ooor oo *•}?■*■»• # * ** -J-J DZh #•»■•«■■»■ * 4J--*4* ■«•■»■■«•-«■ ■«•*** ZC^ZHK ■*-«■*•«• s •!<•**■* 2 2m«Zm **•«••»■ < **■»•# ■«■■»«■■«•■»«• O ■>»»■.":•■»■ I»-|- CXU'J- ■>f •»■ * -K- os #•>;-•*-* •—2" (_)>•-■ ^ ■«•* # 0- * * * » aruj> 3 * ** ■» **■«•* uJCh-oru.02 •«■•»••»*•■»■ a * # < • «• l-ll-3UCZi- * 4** # z ■>>#** <.mO-UJ m»— ■«■•»■■»■ * i— i . ##*4* -J-J21Xf--jZ ***■«■ j ♦ * if * n> #■ * _l —22 *■»•■»■-«■ 3 •k 4^•■H■■K• f-so.rjoo •«•«•■«■-»■ X> fc 4* # 4V »*«« *- UJ ox #**•»• < * * -> ■>> >-«/iZ *-u_ * ** * co ■»- # 4* 4J- OiiOCuu C # * ■»• * ■*• 4J- ■«• * u a: arcr: *■»■•«•■» >H *#■»■* ^Ull/jU-O • •K- ■«■-«■ # w **«•■»■ DI-m XQLL ■»••»■*•* o: * * <• * 2 UJ_J ■«■**■«■ u **■«■* l/"OZ IU.7Q. •* 4<--if * c— *«■*■«• mSO'UU-mih **■»■•«• 2 * ♦ * * r-OOSTO ■*■*-«■ # O *- 41-4*- # 5ZZ!-< Ct2T #*■«■•«• S •ft-*** 0 4f4f ■»••* CU * * ** c uj*-Mm x» r- < ro ► •> ^-"^* «» oo o o.r- o t/>Kl •4- z** mm* »-«3 > 3< < K a v« < •■k •- •> o~ — » r-o o con o «—fn ^f— a:— wO >^ i— 'M >o • QU <*. -> ■— c> O o— OO -~< OC0 ►r- r-or-i — -o — O"^ — ta:oo mc/>0_Jt— Zh> -CL<1-— •t-. ■ — ■ -•CUJ >cm200XhOoi«rcc< Q233X<:0<»( 02T X ooi/)co c/)coco< uj»— mo ZZZZZZZZZZZZUJZ II o UJUJ UJ UJ UL' UJ UJUJ UJ UJ LU LU»-i mm X 00 szz5:z2:z2:sz5;2oa.a< ___M_*_l_MMMM_tOZ a. s: o QQQQOQOQOQOO_Jmm.mQ UJ > o UJ cc <. v X CXLU z UJ X CO r-