Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
The impact of Lagrangian environmental variability on the growth of phytoplankton
(USC Thesis Other)
The impact of Lagrangian environmental variability on the growth of phytoplankton
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
The impact of Lagrangian environmental variability on the growth
of phytoplankton
by
Jessica Zaiss
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(GEOLOGICAL SCIENCES)
December 2024
Copyright 2024 Jessica Zaiss
To all the little girls, from small towns where no one leaves, with an endless curiosity
for how the world works, who dare to dream of becoming a NASA scientist - You can
do hard things.
ii
Acknowledgements
I would like to first and foremost thank my amazing advisor, Dr. Naomi Levine. Thank you for
taking a chance on me, thank you for believing in me when I didn’t believe in myself, thank you
for supporting me through one of the hardest periods of my life, and not giving up on me. Thank
you for guiding me and shaping me into the rigorous scientist that I am today. Thank you, for
making my dreams come true.
I would like to thank my committee members without whom this would not be possible.
Through their guidance, I have confidence that my work is sound, interesting, and thorough. They
have been nothing but supportive and I am grateful for them.
I would like to thank all of my friends. The many, many of them who have been cheering me
on this long journey. Those who have been in my corner since elementary (Bethany!), college,
grad school the first time, grad school the second time, and this third and final time. I could not
have stayed sane (though this is arguable) without my triathlon family; Especially Taryn, Lynne,
and Lauren. Hillary Biscay for her endless support and belief in me. To all my Smashfest Queen
training partners who became friends and then friends who became family.
To Heidi and Rachel - I am so grateful that I did not have to go through losing Mary alone.
You provided love, support, and even just the days of knowing I wasn’t alone in my pain were
appreciated. I am so proud of all of us.
To Nati and Jaime - The two best friends anyone could ever ask for. Thank you for pushing
me on all the days I thought I couldn’t do it; for always believing in me, particularly when I didn’t
believe in myself. Thank you for reminding me to put ”one paragraph in front of the other” and for
providing all the many laughs along the way.
I would like to thank my family. For always supporting me and never asking me when I
was going to get a ”real” job. Also, Rob’s mom and dad, Marcia and Mike, thank you for your
enthusiastic love and support during this adventure.
iii
Lastly, to Rob. For everything. For taking care of the house and dogs during my long work
days. For staying by my side when I was at my worst dealing with the loss of Mary. For always
being supportive of pursuing my dreams.
iv
Table of Contents
Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii
Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii
List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xxi
Chapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
Chapter 2: Impact of Lagrangian Sea Surface Temperature Variability on Southern Ocean
Phytoplankton Community Growth Rates . . . . . . . . . . . . . . . . . . . . . 5
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2.1 Southern Ocean drifter profiles . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2.2 Remote sensing SST . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.2.3 Idealized SST trajectories . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.2.4 Phenotype-based Ecosystem Model . . . . . . . . . . . . . . . . . . . . . 15
2.2.5 Acclimation Rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.3.1 SST variability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.3.2 Impact of variable SSTs on community growth rates . . . . . . . . . . . . 20
2.3.3 Memory Effect . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.3.4 Southern Ocean Drifter Trajectories . . . . . . . . . . . . . . . . . . . . . 23
2.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.4.1 Impact of Acclimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.4.2 Implications for in situ community composition . . . . . . . . . . . . . . . 28
2.4.3 Implications for simulating community growth rates in global biogeochemical models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
2.6 Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
Chapter 3: Lagrangian mixing variablity . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.2.1 Physical Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.2.2 Lagrangian Trajectories . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.2.3 Identification of Submesoscale features . . . . . . . . . . . . . . . . . . . 42
3.2.4 Spatial Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
v
3.3.1 Eulerian Temperature and Mixing Variability . . . . . . . . . . . . . . . . 45
3.3.2 Eulerian Submesoscale Feature Characteristics . . . . . . . . . . . . . . . 46
3.3.2.1 Submesoscale Feature Spatial Patterns . . . . . . . . . . . . . . 46
3.3.2.2 Submesoscale Feature Size . . . . . . . . . . . . . . . . . . . . 49
3.3.2.3 Submesoscale Velocities . . . . . . . . . . . . . . . . . . . . . . 49
3.3.3 Lagrangian Trajectory Dynamics . . . . . . . . . . . . . . . . . . . . . . . 50
3.3.4 Reference Frame Impact on Temperatures Experienced . . . . . . . . . . . 53
3.3.5 Reference Frame Impact on Mixing Rates Encountered . . . . . . . . . . . 53
3.3.6 Reference Frame Impact on Velocity . . . . . . . . . . . . . . . . . . . . . 55
3.3.7 Lagrangian Encounters by Submesoscale Feature Type . . . . . . . . . . . 58
3.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
3.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
Chapter 4: A Lagrangian Lens on Biophysical Interactions in the California Current System 63
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
4.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
4.2.1 Physical Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
4.2.2 Lagrangian Trajectories . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
4.2.3 Phytoplankton Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
4.2.3.1 Biomass and Nitrogen Equations . . . . . . . . . . . . . . . . . 68
4.2.3.2 Cellular C Quota . . . . . . . . . . . . . . . . . . . . . . . . . . 70
4.2.3.3 Nitrogen Uptake Rate . . . . . . . . . . . . . . . . . . . . . . . 70
4.2.3.4 Photosynthesis and Respiration . . . . . . . . . . . . . . . . . . 72
4.2.3.5 Chlorophyll Synthesis and Degradation . . . . . . . . . . . . . . 73
4.2.4 1-D water column model . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
4.2.5 Post-processing model analyses . . . . . . . . . . . . . . . . . . . . . . . 75
4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
4.3.1 Model Validation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
4.3.1.1 Nitrate Concentration . . . . . . . . . . . . . . . . . . . . . . . 76
4.3.1.2 Light, Photosynthesis, and Growth . . . . . . . . . . . . . . . . 79
4.3.2 Biological and Physical Control on Biomass Accumulation . . . . . . . . . 83
4.3.2.1 Biological Controls . . . . . . . . . . . . . . . . . . . . . . . . . 83
4.3.2.2 Physical Controls . . . . . . . . . . . . . . . . . . . . . . . . . . 87
4.3.2.3 Critical Diffusivity . . . . . . . . . . . . . . . . . . . . . . . . . 87
4.3.3 The Impact of Lagrangian Variability . . . . . . . . . . . . . . . . . . . . 89
4.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
4.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
Chapter 5: Conclusions and Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
5.1 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95
5.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
A1 Chapter 1 Supplement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
vi
A1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109
A1.2 Impact of final temperature in the idealized simulations . . . . . . . . . . . 109
A1.3 Impact of reaction norm width in idealized simulations . . . . . . . . . . . 110
A1.4 Model sensitivity to minimum biomass parameter . . . . . . . . . . . . . . 112
A1.5 Comparison of Ecosystem Model Choice . . . . . . . . . . . . . . . . . . 113
A1.5.1 Linear Mortality . . . . . . . . . . . . . . . . . . . . . . . . . . 113
A1.5.2 Quadratic Mortality (used in the main text) . . . . . . . . . . . . 113
A1.5.3 Simple Ecosystem . . . . . . . . . . . . . . . . . . . . . . . . . 114
A1.5.4 Constant grazing . . . . . . . . . . . . . . . . . . . . . . . . . . 115
A1.6 Statistics Calculations for Sensitivity Tests . . . . . . . . . . . . . . . . . . 115
A1.7 Nitrate limitation model . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
A2 Chapter 2 Supplement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148
A2.1 Results sensitivity to initial particle depth . . . . . . . . . . . . . . . . . . 148
A2.2 Vertical Velocity Filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148
A2.3 Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150
A3 Chapter 3 Supplement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156
A3.1 Figures and Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156
vii
List of Tables
A1.1 Results of SSTmax variability analysis . . . . . . . . . . . . . . . . . . . . . . . . 118
A1.2 Estimates of acclimation times (lower bounds) for the polar diatom F. cylindrus
for lab-cultures initiated at 2.5 °C and then exposed to first to either cooling or
warming of ≈1.5 °C and 3 °C for the first acclimation, and then subsequent step
increases or decreases in temperature, as part of the preparation for a thermal
reaction norm experiment. The estimates growth rates were used in conjunction
with the change in temperature to compute the range of acclimation times
presented. For further information contact philip.boyd@utas.edu.au. . . . . . . . . 119
A3.1 Biology Model Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157
A3.2 Physical Model Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157
A3.3 Initial Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158
viii
List of Figures
2.1 Lagrangian versus Eulerian reference frames. Lagrangrian reference frames
follow the water parcel itself through time. The Eulerian reference frame refers to
a fixed point in space (e.g. buoys or mooring stations) where advection of water
parcels floating past the fixed point generates temporal variability. Panels (a), (b),
and (c) depict three different water masses (grey, blue, and green) as they each
pass through the fixed Eulerian location (red dot) at times t1, t2, and t3¬. Panel
(d) shows the temperature of each water mass through time (grey, blue, and green
lines) as well as the temperature recorded at the Eulerian location (red line). Note
in panel (d) that overall temperature variability in the Lagrangian reference frame
(grey, blue and green lines) is much greater than that in the Eulerian reference
frame (red dots) though this may not always be the case. . . . . . . . . . . . . . . 10
2.2 The impact of SST variability on individual phenotype growth rate. (a) The
temperature related growth response for a phenotype with a skewed shaped
reaction norm. The values for the optimum growth temperature (Topt) and the
corresponding maximum growth rate (µmax) are shown with dashed lines. (b) The
90-day SST profile of an example drifter trajectory (black) and the associated
changes in phenotype growth rate (blue). The orange and red arrows in the
top panel indicate the change in the phenotype growth rate associated with the
corresponding changes in SST in the bottom panel. . . . . . . . . . . . . . . . . . 11
2.3 a) Map of all 90-day drifter trajectories (n = 2190) colored by SST. Two example
trajectories are highlighted in purple and magenta. b) Reaction norms for each
of the 319 phenotypes in the ecosystem model. The grey lines represent all the
phenotype reaction norms and the green lines are example phenotypes to highlight
the reaction norm shape. c and d) Example trajectories and their resulting model
outputs. The top panels show the SST (colors), the community growth rate
estimated using the Eppley curve (dashed line), and the community growth rate
from our phenotype-based model as calculated using Eq. 5 (solid line). The
bottom panel shows the biomass through time of each phenotype (grey lines). The
blue line follows the phenotype with the highest initial biomass, the red dashed
line follows the phenotype that has the highest biomass at the end of the 90 days,
and the green line follows the phenotype that has a Topt equal to the mean SST of
the trajectory. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.4 SST variability analysis. The frequency of ∆SSTmax changes from the drifter
segments over different window lengths are shown. Data are presented as total
percent of data that fall within that window length such that each row sums to
100%. There is a general pattern of increasing magnitudes of ∆SSTmax over longer
window lengths. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
ix
2.5 Simulated response of phytoplankton community with skewed shaped reaction
norm to increasing ∆SST (see Appendix for decreasing ∆SST conditions and
broad reaction norm results). Panel (a) plots the decline in community growth rate
in the phenotype model that results from the SST moving out of the thermal niche
of the original population (see Methods and Figure A1.16). Data that are greyed
out represent ∆SST and window length combinations that were not supported
by the results from Figure 3. Panel (b) shows the percent difference between
the Eppley growth model approximation and the phenotype modeled community
growth rates at the point where SST stabilizes (see Figure A1.16 for example).
Panel (c) plots the memory effect length associated with SST changes in the
idealized simulations. This represents the time it takes for the community growth
rate to be within 5% of the steady state growth rate at the final SST from the first
time-step that SST is constant (See Figure A1.16 for example). . . . . . . . . . . 26
2.6 Impact of SST variability on community growth rate. The average percent
difference in community growth rate between the phenotype model and the Eppley
growth model from the 90-day drifter segments are plotted against the standard
deviation (1σ) of the drifter SST. Each segment is colored by the mean SST.
Results from the idealized trajectories are shown as black diamonds with filled
circles denoting increasing SST trajectories and open circles denoting decreasing
SSTs. Pink triangles represent the two example trajectories from Figure 3. Results
shown here are for skewed shaped reaction norms, see Figure A1.20 for results for
the broad shaped reaction norms. . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.7 The impact of acclimation on the number of generations for which the effect
persists (memory length). Acclimation rates that were slower than the rate of
SST change resulted in longer memory lengths than for simulations in which
acclimation rate was equal to or faster than the SST rate of change. . . . . . . . . 29
2.8 The impact of Lagrangian and Eulerian variability on community composition.
Here we plot the difference between the Topt of the most abundant phenotype
at the end of each 90-day trajectory and the final SST for the drifter trajectory
(x-axis) and the satellite data (y-axis). The final SSTs for the drifter and satellite
data are not statistically different (t-test, 95% CI). Therefore, deviations from the
1:1 line demonstrate the impact of a Lagrangian versus Eulerian reference frame
on community composition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
x
2.9 Distribution of SST variability (a-c) and the deviation in community growth rate
from the Eppley growth model (d-f) over the Southern Ocean (>30 °S). Only
those drifters which overlap in space and time with the satellite data are shown.
For full results, see Figure A1.29. Three key regions of high SST variability
stand out: Malvinas-Brazil confluence zone, the Agulhas Retroflection, and the
Subtropical front. These regions have enhanced SST variability in all datasets but
higher variability in the Lagrangian trajectories. These high variability regions
correspond to large differences between the phenotype model growth rates and
the Eppley approximation of growth, a pattern consistent across all three sets of
simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.1 Examples of Lagrangian particle trajectories. a) Starting locations of the Lagrangian
particles at the first timestep, colored according to starting latitude. The background is
colored by sea surface temperature [] at the first time step. b) Final location of the particles
on day 66. The color represents the latitude that the particles started at as shown in (a).
The size of the particle represents the depth of the particle with larger sizes indicating
deeper depths as specified in the legend (meters). c) Trajectory of an example Lagrangian
particle through time. The background is colored by the relative vorticity at day 54.5
and the trajectory is colored by simulation day. At simulation day 54.5, the Lagrangian
particle encounters a cyclonic frontal feature, marked by a star. (d) The depth through
time of the example trajectory shown in (c). The line is colored by the relative vorticity
that the Lagrangian particle encountered. The initial encounter with the frontal feature is
marked by a start at simulation day 54.5 after which the particle is rapidly subducted by
the frontal flow. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.2 Model temperature and mixing conditions over the simulation period. a) Mean temperature
at 39.72 5m depth for the period of March 27th -June 1st . b) standard deviation (2σ)
at 39.725 m depth. c) Mean diffusivity at 39.725 m for the whole simulation period. d)
diffusivity standard deviation (2σ) divided by the mean diffusivity. . . . . . . . . . . . . 47
3.3 Prevalence of submesoscale features. a) As a proportion of the surface area of each region
b) As a proportion of Lagrangian particle encounters per region. . . . . . . . . . . . . . 48
3.4 Mean vertical and horizontal velocities as calculated from the Eulerian reference frame
compared the mean velocities encountered in the Lagrangian reference frame. The solid
black line indicates the Eulerian mean equal to the Lagrangian mean. . . . . . . . . . . . 50
3.5 Lagrangian particle trajectories over the entire simulation period. Trajectories are colored
by release latitude (see Figure 1 for release locations and final depths). . . . . . . . . . . 51
3.6 Mean mixed layer depth across all Lagrangian particles over time (blue). Blue shading
represents the range of the mixed layer depth across all particles at each timestep. Also
shown is the mean depth of the Lagrangian particles over time (black) and the maximum
depth across all particles at each time step (dotted black line). The green line represents
the proportion of particles that are within the mixed layer. Changes in the proportion
of particles in the mixed layer results from mixed layer dynamics opposed to vertical
advection of the particles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
xi
3.7 Difference between Lagrangian and Eulerian reference frames as a function of particle
advection. (a) As Lagrangian particles are advected offshore, the mean temperature
along the trajectory decreases relative to the Eulerian mean. However, the temperature
variability experienced by the Lagrangian particles advected offshore is higher than the
corresponding Eulerian variability. (b) Same as panel (a) but showing differences in
diffusive mixing values between the two reference frames. Although the difference in
means between the two reference frames follows the same spatial pattern as temperature,
the Lagrangian variability decreases relative to the Eulerian reference frame with offshore
advection. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.8 Kernel density estimates (analogous to probability distributions) of the time spent in a
submesoscale feature by the Lagrangian particles. The data are separated by region and
data type. Generally, particles spent the longest time in cyclonic features. However, the
mean time spent in a feature is 1.6 hours, much shorter than the lifetime of the features. . . 56
3.9 Comparison of Eulerian occurrence rate by area to the Lagrangian encounter rate by time
spent in the mixed layer. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.10 Proportion of total vertical advection resulting solely from advection due to submesoscale
features. The analysis was performed across all Lagrangian particles. Lagrangian
trajectories with submesoscale proportions of less than 1.9*10−6
are shown. . . . . . . . 60
4.1 : Examples of Lagrangian drifter trajectories. A) Starting locations of all the Lagrangian
drifters at the first time step. drifters were initialized at a depth of 100m and are colored
according to starting latitude. The background is colored by sea surface temperature [] at
the first time step. B) Final location of the drifters on day 66. The color represents the
latitude that the drifters started at in (A). The size of the drifter represents the depth of
the drifter. Larger point size mean deeper drifters as specified in the legend (meters). C)
Trajectory of a Lagrangian drifter through time. The background is colored by the relative
vorticity at day 54.5 and the trajectory is colored by simulation day. At simulation day
54.5, the Lagrangian drifter encounters a cyclonic frontal feature, marked by a star. (D)
The depth through time of the trajectory in (C). The line is colored by the relative vorticity
that the Lagrangian drifter encounters. The initial encounter with the frontal feature is
marked by a start at simulation day 54.5. The drifter is rapidly subducted by the frontal flow. 67
4.2 : Example of Lagrangian phytoplankton model. (Top) Idealized latitude and longitude
of a Lagrangian particle over the whole simulation. (Middle) Depth profile through time
of the idealized trajectory in the top panel. (Bottom) 1-D water column model with
phytoplankton biology advected along the Lagrangian particle trajectory over the entire
simulation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
4.3 The relationship between nitrate and depth for the CalCOFI bottle data (light blue), our
model data (dark blue), and measurements from Li et al. (2010) (orange diamonds). The
yellow and blue lines represent the bottle and model means, respectively. Dashed lines of
the same colors represent the 2σ standard deviation. . . . . . . . . . . . . . . . . . . . 77
xii
4.4 The relationship between nitrate and temperature for the bottle date, the Li et al., 2010
measurements, our model data, and the Lucas model (Lucas et al., 2011). Solid dark blue
points and light blue points are the model and CalCOFI bottle data, respectively. The
yellow and blue lines represent the bottle and model means, respectively. Dashed lines
of the same colors represent the 2σ standard deviation. Orange diamonds represent the
measurements from Li et al., 2010. The gray solid line shows the temperature-nitrate
relationship Lucas et al., 2011b. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
4.5 Vertical profile of mean PAR in the model compared to observed PAR data from Mitchell
et al., 2017. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
4.6 Vertical profiles of Chl with depth. Data from CalCOFI (light blue dots, yellow lines), our
model (dark blue dots, green lines), and Li et al. (2010) (orange diamonds) are shown. . . 80
4.7 Realized growth rates from the model (dark blue dots) compared to estimated growth rates
from Li et al. (2010) (orange dots) as a function of depth. The model mean( green line)
and standard deviation (dotted lines, 2σ) encompass 95% of the Li et al. (2010) data. . . . 81
4.8 Vertical profiles of carbon from the model (dark blue dots) and Li et al. (2010). The model
captures both the extreme and general profile of carbon from the Li et al. data. . . . . . . 82
4.9 Proportion of time that changes in biomass are controlled by biological and physical
processes. Data are ordered from highest proportion of biology dominated changes to
lowest. Biological controls on biomass change are light limitation, nutrient limitation,
and temperature limitation. When biological processes control the change in biomass, the
majority of that time nutrient limitation is the dominant factor. Physical processes that
impact biomass concentration are advection and diffusion. Diffusion is the most dominant
physical process. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
4.10 When biology was the controlling process, the change in biomass was on average positive
(0.18 mmol C/m3
/hr or 4.32 mmol C/m3
/day). When diffusive mixing was the controlling
process, the change in biomass skewed negative with an average change of -0.19 mmol
C/m3
/hr (4.6 mmol C/m3
/day) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
4.11 Limiting biological factors when changes in biomass are due to biological processes. (a)
When biomass concentration increases due to biological processes, nutrients are typically
the least limiting factor, although relieving light limitation is also important. (b) When
biomass concentration decreases due to biological processes, it is approximately half the
time due to light limitation while the other half is due to nutrient limitation. . . . . . . . . 86
4.12 Dominating physical process when the change in biomass is a) positive, b) negative.
Turbulent/diffusive mixing is the controls biomass most often but advection does
sometimes lead to biomass increases. . . . . . . . . . . . . . . . . . . . . . . . . . . 88
4.13 Relationship between diffusivity (Kz) and the change in biomass per time step. Colors
of the data points represent the growth rate. Particles with the highest growth rates are
counter intuitively most associated with biomass decreases. . . . . . . . . . . . . . . . . 89
xiii
4.14 Two example trajectories over the 66.7 day simulation. Stars mark the start of the
simulation. These particles started 395km apart and ended only 3.1km from one another. . 90
4.15 Comparison of two example trajectories. (a,b) Biomass concentration through time on
the left axis and on the right axis are the diffusive mixing rates. Shaded regions represent
timesteps in which the change in biomass was dominated by physical processes. (c,d)
The biomass concentration along the entire 200m water column for the whole simulation
period. Also shown are the depths of the particle along the trajectory, the mixed layer, and
the base of the euphotic zone. (e,f) Same as in (c,d) but nitrate concentration instead of
biomass. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
4.16 Comparison of the whole water column nitrate (a) and biomass (b) concentrations of ET1
and ET2 at the final timestep of the simulation. . . . . . . . . . . . . . . . . . . . . . 92
A1.1 Box plots of the percent of the SST variability in the drifter trajectory that is
accounted for by the smoothed spline. Each of the 2,190 90-day drifter and spline
trajectories was broken up into windows in 1-day increments from 1 to 90 days.
The standard deviation of the drifter trajectory is the sum of the standard deviation
of the smoothed spline plus some noise term. From this, the variability accounted
for by the spline for each window, for each trajectory was recorded with the results
shown. As expected, over longer window lengths the spline accounts for higher
percentage of the overall variability. . . . . . . . . . . . . . . . . . . . . . . . . . 120
A1.2 The impact of final SST on percent difference between the individual based
model and the Eppley growth model, relative to Eppley growth ([Eppley –
phenotype]/Eppley). The results from the simulations in the main text are
compared to simulations with final SSTs of 10 °C (a,c) and 20 °C (b,d) for both
the skewed (top row) and broad (bottom row) shaped reaction norms. Open
data points represent decreasing ∆SSTs and filled in data are increasing ∆SSTs.
The black line indicates the 1-1 line. There is no statistical difference between
simulations with differing final SSTs (95% CI). . . . . . . . . . . . . . . . . . . . 121
A1.3 The impact of final SST on memory length. The results from the simulations in
the main text are compared to simulations with final SSTs of 10 °C (a,c) and 20
°C (b,d) for both the skewed (top row) and broad (bottom row) shaped reaction
norms. Open data points represent decreasing ∆SSTs and filled in data are
increasing ∆SSTs. The black line is the 1-1 line. There is no statistical difference
between simulations with differing final SSTs (95% CI). . . . . . . . . . . . . . . 122
xiv
A1.4 The impact of reaction norm width on percent difference from the Eppley
growth model. The results from the simulations in the main text are compared
to simulations with narrower (a,c) and wider (b,d) reaction norms for both the
skewed (top row) and broad (bottom row) shaped reaction norms. The black line
is the 1-1 line. Closed data points represent increasing ∆SSTs and open circles
represent decreasing ∆SSTs. For broad reaction norms, increasing the reaction
norm width increases the difference between the phenotype model and the Eppley
growth model. There was no significant difference between the simulations for
skewed reaction norms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
A1.5 The impact of reaction norm width on memory length. The results from the
simulations in the main text are compared to simulations with narrower (left
column) and wider (right column) reaction norm widths for both the skewed (top
row) and broad (bottom row) shaped reaction norms. The black line is the 1-1
line. Closed data points represent increasing ∆SSTs and open circles represent
decreasing ∆SSTs. Broad reaction norms are most affected by reaction norm
width, but increasing the reaction norm width had a significant (95% CI) impact
on the memory length for the skewed reaction norms as well. For small ∆SSTs
(2-3o°C), narrower reaction norms have longer memory lengths. When ∆SSTs are
large (8-9 °C), the memory length is shorter for communities with more narrow
reaction norms. For moderate ∆SST changes, (4-6 °C), the width of the reaction
norm has minimal impact on the memory length. The memory lengths associated
with moderate ∆SSTs are typically the longest memory lengths, just as in the main
text. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
A1.6 The impact of minimum biomass on deviation from Q10. The results from the
simulations in the main text (x-axis) are compared to simulations with an order
of magnitude smaller minimum biomass for both skewed (left) and broad (right)
shaped reaction norms. The black line is the 1-1 line. Filled in data points
represent increasing ∆SSTs and open data points are decreasing ∆SSTs. The
minimum biomass impact is significant at the 95% CI with an average increase in
offset from Q10, for both reaction norm shapes. . . . . . . . . . . . . . . . . . . . 125
A1.7 The impact of minimum biomass on memory length. The results from the
simulations in the main text (x-axis) are compared to simulations with an order
of magnitude smaller minimum biomass for ‘both skewed (left) and broad (right)
shaped reaction norms. The black line is the 1-1 line. Filled in data points
represent increasing ∆SSTs and open data points are decreasing ∆SSTs. The
minimum biomass impact is significant at the 95% CI with an average increase in
memory length for both reaction norm shapes. However, the pattern of moderate
∆SSTs exhibiting the longest memory effects were robust across all simulations. . 125
xv
A1.8 Comparison between different ecosystem model results for community growth for
an idealized simulation with an increase of 4 °C over 21 days. For community
growth rates, all models show similar qualitative results indicating a decrease in
growth rate over the transient conditions culminating in a growth rate minimum
when SSTs stabilize. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
A1.9 Nitrate-temperature relationship. Relationship between temperature and nitrate in
the Southern Ocean as measured by Bio-ARGO floats. Data was collected from
the Southern Ocean Carbon and Climate Observations and Modeling (SOCCOM)
Project funded by National Science Foundation, Division of Polar Programs (NSF
PLR -1425989), supplemented by NOAA and NASA. . . . . . . . . . . . . . . . . 127
A1.10Nitrate limitation versus temperature. Due to the non-linear relationship, for
temperatures greater than 14 °C, nitrate concentrations minimally impact growth
rates. Once SST goes above 14 °C, the limitation rapidly decreases to the
minimum which results in large growth rate decreases. . . . . . . . . . . . . . . . 128
A1.11Examples of the impact of nitrate limitation on community growth rates. When
SSTs are below 14 °C, the growth rate with nitrate limitation and that without
limitation are nearly identical (Example Trajectory 1). As growth rates approach
14 °C, the growth rate dynamics become more variable but follow the same pattern
as that with temperature only limitation (Example Trajectory 2). When SSTs
exceed 14 °C, the growth rate dynamics are muted due to lower overall growth
rates but still mimic those of the temperature limited only growth rates (Example
Trajectory 3). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
A1.12Probability density functions of the absolute value of the maximum change in SST
over 7, 21, 45, and 90 days for the drifter trajectories (blue) and the smoothed
splines of the trajectory SSTs (red). . . . . . . . . . . . . . . . . . . . . . . . . . 130
A1.13The standard deviation (1σ) as a function of ∆SSTmax over different ∆tmax
windows. ∆SSTmax drives the variability across the ∆tmax window lengths. . . . . . 131
A1.14Results from SST variability analyses for the drifter (left) and satellite trajectories
(right) showing most common SST changes for each time window. Data are
presented as total percent of data that fall within that ∆SST bin for the window
length. Each row sums to 100%. Although the magnitudes of variability are
similar, the nature of that variability is different with the Lagrangian reference
frame experiencing more variability consistent with longer memory effects. . . . . 132
A1.15Daily rates of SST change for drifter trajectories. The rates of change were
calculated as the range of the recorded SST values over a 1-day moving window
for a total of n = 781,749 data points for 197,100 days. . . . . . . . . . . . . . . . 133
xvi
A1.16Community growth rates for each of the 100 simulations (grey lines) for an
increase of 4 °C over 7, 21, 45, and 90 days for skewed shaped reaction norms.
The black line is the Q10 simulated community growth rate and the blue line is
the SST profile for the simulation. The locations marked a1 and a2 in the 21-day
panel represent the timesteps used to calculate the percent change in growth rate
associated with transient SSTs as shown in Figures 2b. This metric was calculated
as (a1 – a2)*100/a1. The locations marked b1 and b2 in the 45-day panel represent
the timesteps used to calculate the percent difference in growth rates between the
Q10 parameterized growth and the phenotype model as shown in Figure 2c, S11.
This metric was calculated as (b1 – b2)*100/b1. The locations marked c1 and c2
in the 90-day panel point to the timesteps used to calculate the memory length.
The dashed grey lines represent ±5% of the final, stabilized community growth
rate which was used as the threshold for the memory effect which was defined
as the time in days between c1 when SSTs stabilize and c2 when the community
growth rate crosses the threshold. . . . . . . . . . . . . . . . . . . . . . . . . . . 134
A1.17Full results for the change in growth rate in the idealized simulations for the
skewed shaped reaction norms (top row) and the broad shaped reaction norms
(bottom row) under both increasing SSTs (left column) and decreasing SSTs (right
column). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
A1.18Full results for the percent difference from the Eppley growth approximation at
the time step when SSTs stabilize in the idealized simulations for the skewed
shaped reaction norms (top row) and the broad shaped reaction norms (bottom
row) under both increasing SSTs (left column) and decreasing SSTs (right column). 136
A1.19Full results for the length of the memory effect in the idealized simulations for
the skewed shaped reaction norms (top row) and the broad shaped reaction norms
(bottom row) under both increasing SSTs (left column) and decreasing SSTs (right
column). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
A1.20The 90-day average percent difference between community growth rates
determined via the Q10 method and the phenotype-based model versus standard
deviation (1σ) of SST over the 90-day trajectory. Drifter data are represented by
circles colored according to their mean SST. Black diamonds represent the first
90 days of the idealized trajectories; filled diamonds are the idealized trajectories
for which SSTs increase and open black diamonds are idealized trajectories with
decreasing SSTs. Pink triangles represent the two example trajectories from
Figure 1 in the main text. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
A1.21Comparison of mean community growth rate over the entire 90-day trajectory for
the real trajectories and their spline simulations for skewed (left) and broad (right)
shaped reaction norms. With each reaction norm shape, smoothing the small-scale
noise did not impact overall biomass-weighted community growth rates (95% CI,
t-test) further supporting that small-scale noise does not induce a memory effect. . 139
xvii
A1.22Example of impact of acclimation. (Top) Example idealized SST trajectory of
changing 2 °C in 7 days with acclimation rates of 0.2 °C day-1 and 0.3 °C day-1.
Other acclimation rates not shown as they plot along the No Acclimation line
because those rates are faster than the rate of change. (Middle) Community
growth rates for skewed reaction norms over each of the simulations for the no
acclimation simulations (grey lines) and the simulations with an acclimation
rate of 0.2 °C day-1 (green lines). Dashed lines represent the thresholds used to
calculate the memory length. (Bottom) Same as the middle panel but for broad
reaction norms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
A1.23Impact of acclimation on memory length on the skewed reaction norms (top row)
and the broad reaction norms (bottom row) in both the increasing (left column) and
decreasing (right column) ∆SST directions. Dashed lines represent the memory
lengths calculated for the simulations that did not incorporate acclimation. When
acclimation rates are greater than or equal to the rate of SST change, there is no
difference between the simulations that incorporated acclimation and those that
did not. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
A1.24Comparison of the final SST for the drifter and the satellite data. The data from
both sources represent the same location in space and time so the data should be
similar and in fact, are not statistically different from one another (ttest, 95% CI).
The grey line represents the 1-1 line. . . . . . . . . . . . . . . . . . . . . . . . . . 142
A1.25Same as Figure 4 in the main text but for broad shaped reaction norms. The
impact of Lagrangian and Eulerian variability on community composition. Here
we plot the difference between the Topt of the most abundant phenotype at the
end of each 90-day trajectory and the final SST for the drifter trajectory (x-axis)
and the satellite data (y-axis). The final SSTs for the drifter and satellite data
are statistically identical (t-test, 95% CI). Therefore, deviations from the 1:1
line demonstrate the impact of a Lagrangian versus Eulerian reference frame on
community composition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143
A1.26The impact of SST variability on community composition. (Top) An example
90-day drifter trajectory and the satellite SST data for the final location of the
drifter over the same 90 days shown as solid lines. The dashed lines are the Topt
of the most abundant phenotype at each timestep. (Bottom) The biomass of each
phenotype with a skewed shaped reaction norms at day 90 for the satellite and
drifter trajectories. In this example, the offset between the final SST is -0.60 °C
for the drifter and -0.09 °C for the satellite data. The difference in the magnitude
of the offset between the two data sets represents the difference in the variability
of the SSTs. However, in this example, the satellite SSTs stay relatively constant
whereas the drifter SSTs experience a rapid increase of 3.5 °C in 4 days beginning
Sept. 29. Because the drifter SSTs remain relatively constant through the end of
the 90 days, the community is able to adjust to the new environment before the
end of the simulation which results a community Topt that reflects the SSTs at day
90 for both the satellite and the drifters. . . . . . . . . . . . . . . . . . . . . . . . 144
xviii
A1.27The impact of SST variability on community composition. (Top) An example
90-day drifter trajectory and the satellite SST data for the final location of the
drifter over the same 90 days shown as solid lines. The dashed lines are the Topt
of the most abundant phenotype at each timestep. (Bottom) The biomass of each
phenotype with a skewed shaped reaction norms at day 90 for the satellite and
drifter trajectories. In this example, the offset between the final SST is -0.23 °C
for the drifter and -3.1 °C for the satellite data. The difference in the magnitude of
the offset between the two data sets represents the difference in the variability of
the SSTs. Here, the drifter SSTs gradually increase over the 90 days which allows
the community to continuously track the changes in SST whereas the satellite
SSTs are relatively stable and then rapidly decrease from 17.7 °C on March 10 to
13.8 °C on March 17. Due to the long memory effect associated with this rate and
magnitude of change, the community was not able to track the SST change which
resulted in a large offset between the final SST and the Topt of the most abundant
phenotype at day 90. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
A1.28Percent reduction in community growth rate from the phenotype model from the
Eppley growth model for the drifters, the satellite derived trajectories, and the
satellite point data. In the Lagrangian reference frame (both sets of trajectories)
the community growth rate from the phenotype model is lower than in the Eulerian
reference frame (satellite point data) as a result of the different SST variability
encountered in each reference frame. . . . . . . . . . . . . . . . . . . . . . . . . . 146
A1.29Full drifter results from Figure 9 in the main text. This figure includes drifter
trajectories that do not overlap in time with the satellite data. The spatial patterns
are the same as in the main text. . . . . . . . . . . . . . . . . . . . . . . . . . . . 147
A2.1 Spatial Occurrence rate of submesoscale features versus particle encounter rate for
Lagrangian particles initialized at 100m, 50m, and 10m. Closed data points represent data
from the onshore regions. Open data points are from the offshore regions. Triangles are
for the Southern regions and squares represent the northern regions. . . . . . . . . . . . 148
A2.2 The impact of filtering out large-scale vertical velocities to isolate submesoscale features.
(left) Raw, unfiltered vertical velocity at a depth of 9.67m. (middle) Filtered vertical
velocities at 9.67m depth. Submesoscale filaments and eddies are revealed. (right)
Rossby number at 9.67m. Submesoscale features are identified by absolute values >0.75.
The locations of submesoscale features aligns with the vertical velocities revealed after
filtering the raw data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150
A2.3 Standard deviation (2σ) of diffusivity over the whole simulation for a depth of 39.725 m. . 151
xix
A2.4 The size of submesoscale features varied significantly among regions and feature types.
Cyclonic eddies were consistently the largest, averaging 2.36 ±0.73 km2
(2σ), followed
by fronts (1.32 ±0.43 km2
), and anticyclonic eddies (1.15 ±0.42 km2
) (Figure S2). The
smaller size for anticyclonic eddies aligns with regional observations (Payandeh et al.,
2023) and is likely attributed to their susceptibility to instabilities that limits their growth
(Shi et al., 2024). Notably, features in the southern regions were, on average, 29% larger
than their northern region counter parts. Despite some overlap in size distributions, a
two-sided t-test confirmed significant differences in mean feature size among regions and
feature types (p<0.05). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152
A2.5 Distribution of submesoscale feature and background vertical velocities. (Left) Eulerian
distribution (Right) Lagrangian distribution. The mean vertical velocities were statistically
different form each other (two-sided t-test, p<0.05) and from the background average. . . 153
A2.6 Distribution of submesoscale feature and background horizontal velocities(Left) Eulerian
distribution (Right) Lagrangian distribution . . . . . . . . . . . . . . . . . . . . . . . 153
A2.7 Distribution of mean Eulerian vertical velocities for submesoscale features and the
background. Data represent the filtered vertical velocities to remove large-scale
contributions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154
A2.8 Distribution of Lagrangian vertical velocities for submesoscale features and the
background. Data represent the filtered vertical velocities to remove large-scale
contributions. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155
A3.1 Schematic of the quota model and the 1D, Lagrangian water column model . . . . . . . . 156
A3.2 Impact of different initial conditions. Two different sets of initial conditions were
compared to the original initial conditions used in the main text. This sensitivity test was
performed on all 295 particles released at 38 . This histogram shows the days that it took
for the biomass of the two sensitivity test conditions to equal that of the original initial
conditions. From this, we determined that initial conditions largely only impact the first
10 days of the trajectories and so those are excluded from analyses. . . . . . . . . . . . . 158
xx
Abstract
Phytoplankton are vital to the global carbon cycle, contributing half of the world’s photosynthesis. This process removes carbon dioxide from the atmosphere, helping regulate Earth’s climate.
However, large-scale climate models struggle to represent phytoplankton dynamics due to their
coarse spatial resolution ( 10-150km) which smooths spatial variability. Smoothing spatial variability overlooks the nonlinear response of phytoplankton growth and biomass accumulation rates
to environmental changes. This dissertation investigates the critical bio-physical interactions that
occur in the ocean that govern the growth of phytoplankton and resulting biomass accumulation.
I focus on the timescales ranging from hourly to seasonal dynamics and spatial scales ranging
from 1km to the entire Southern Ocean. To accurately determine the bio-physical response to
environmental fluctuations, it is essential to adopt the phytoplankton’s reference frame. As phytoplankton are largely passive drifters, they experience environmental variability in the “along-flow”
reference frame (Lagrangian reference frame), which is distinctly different from the fixed-point,
Eulerian reference frame. Using a combination of Lagrangian based methods and biogeochemical
models, I quantify how phytoplankton growth and biomass accumulation respond to varying rates
and magnitudes of temperature changes over hours to seasons. I also unravel the complex interactions between advection/diffusion and biomass accumulation at small scales (hourly, 1km). I
show that moderate temperature changes (3-6 oC) within 7-21 days cause lasting decreases in community growth rates, a phenomenon not captured by typical, large-scale model, growth dynamics.
Furthermore, resolving fine-scale physics reveals that biomass accumulation is driven by a complex interplay between mixing and non-linear growth dynamics. These findings challenge existing
hypotheses on biomass accumulation and highlight the intricate dynamics governing phytoplankton populations. This improved understanding provides a basis for refining large-scale models,
ultimately enabling more accurate predictions of the ocean’s role in the global carbon cycle and
climate regulation.
xxi
Chapter 1
Introduction
Ocean phytoplankton are a crucial part of the global carbon cycle. They contribute 50% of global
photosynthesis (Field et al., 1998) and are the beginning of the biological pump which acts to
remove carbon dioxide from the atmosphere and sequester it in the deep ocean for up to a millennia
(Falkowski et al., 2003). The inclusion of phytoplankton in global climate models is imperative
for studying how our climate will respond to ever-increasing atmospheric carbon dioxide.
Over the past few decades, great strides have been made in representing phytoplankton dynamics in large-scale global models. Most global climate models now include at least one type
of phytoplankton (Sef´ erian et al., 2020). However, global climate models have enormous compu- ´
tation demands which limits the spatial resolution at which they can be run. This has important
implications for the ability of these models to accurately capture phytoplankton dynamics. Physical processes that impact phytoplankton growth and biomass accumulation dynamics vary from
large, mesoscale eddies ( 100km, months - years) to small, submesoscale fronts ( 1-10km, hoursdays) Landry et al., 2009; Li et al., 2012; Mangolte et al., 2022; Gaube et al., 2013; Chelton et al.,
2011. Spatial patterns of phytoplankton accumulation reflect the cumulative effect of large- and
small-scale environmental variability such as variations in temperature, nutrients, vertical and horizontal advection, and diffusive mixing rates. Within models, everything at the sub-grid cell level
(i.e. all model state variables such as temperature and nutrients) are assumed to be homogenized,
mean values.Thus, for global climate models, phytoplankton experience environmental conditions
that are averaged over tens of square kilometers. In addition to not resolving fine-scale spatial gradients, this also reduces the temporal variability experienced by plankton. By nature of non-linear
growth dynamics, the growth simulated in an average environmental condition will not be the same
1
as the mean growth rate under variable conditions (Jensen, 1906). Better representation of phytoplankton dynamics in coarse-grid simulations requires a better understanding of phytoplankton
dynamics at fine spatial and temporal resolutions.
While the impact of large, mesoscale dynamics on phytoplankton growth and accumulation
have been well studied, further research is needed to elucidate the biophysical interactions at the
submesoscale. Mesoscale physical processes are dominated by largely horizontal, advective flow,
a consequence of the geostrophic balance between the horizontal pressure gradient force and the
Coriolis effect. As such, these large-scale features are stable and can persist for many months, up to
years. Submesoscale dynamics represent balanced flows that are susceptible to disturbances. The
fine-scale fronts and eddies often form because of mesoscale strain between two large-scale eddies.
In the region between two mesoscale eddies, the strain from the large-scale dynamics often creates
an area of convergence in which horizontal buoyancy gradients form. As the strain field persists
and the buoyancy gradient increases, the delicately balanced density gradient becomes unstable
and rapidly adjusts to generate more stratified conditions. Specifically, water on the light side of
the gradient is quickly upwelled while the dense side subducts. The vertical velocities associated
with these features can reach as high as 100 m/day (Mahadevan et al., 2006).
The rapid vertical advection associated with submesoscale features has the potential to significantly impact microbes living in the surface ocean. The rapid upwelling and subsequent water column stabilization injects nutrients and light limited phytoplankton into the euphotic zone.
Once there, the stratification acts to retain them in the sunlit zone enabling growth. Conversely,
downwelling removes biomass from the surface waters. However, it has also been suggested that
because the rapid vertical velocities associated with submesoscale features are so large that biological impacts may be negligible due to rapid vertical exchange into and out of the euphotic zone
(Levy et al., 2012a). Furthermore, modeling studies demonstrate that submesoscale features have ´
the largest impact when they extend past the mixed layer (Levy et al., 2018). Thus, while remote ´
2
sensing, observations, and high-resolution model simulations show phytoplankton spatial heterogeneity associated with these fine-scale physics, the relative importance and time-scales associated
with these different processes is still unclear.
Extracting the impact of spatial and temporal variability on phytoplankton growth dynamics is
further complicated by the reference frame in which phytoplankton are simulated and observations
are collected. Phytoplankton are largely passive drifters and, as such, experience the environment
along the flow of the current, or in the Lagrangian reference frame. Conversely, model simulations
and most data collection methods are in the fixed-point, Eulerian reference frame (e.g. stationary
time series measurements). The mean environmental conditions, and environmental variability,
have been shown to differ between the two reference frames (Doblin et al., 2016; Kida et al., 2017;
Kling et al., 2019).
This dissertation addresses critical knowledge gaps in the time-scales of phytoplankton responses to fine-scale environmental variability by investigating the bio-physical interactions governing phytoplankton growth and biomass accumulation from a Lagrangian perspective, with a
particular focus on fine-scale physical processes.
In my first chapter, I quantified the impact of temperature variability on phytoplankton community growth and biomass accumulation dynamics. I used a suite of idealized and realistic simulations on timescales ranging from hourly to seasonal and employed both Lagrangian and Eulerian
methods to quantify the impact of reference frame choice on overall phytoplankton growth rates.
In my second chapter, I looked further into the physical processes that can cause spatial and
temporal variability. I evaluated how often a passively drifting, Lagrangian particle encounters
small-scale, submesoscale features, and how long the particle stays within the features. I also
compared mean environmental conditions between the Lagrangian and Eulerian reference frames
For my third chapter, I developed a novel, Lagrangian water column model that simulates the
upper 200m of the water column along a Lagrangian trajectory. I used this model to evaluate
the dynamics between biomass changes due to biology (e.g. growth and mortality) and physical
processes (e.g. advection and mixing). For this work, I utilized output from a high resolution
3
physical model (500m horizontal resolution, hourly output) (Torres et al., 2022). I analyzed the
importance of limiting factors in biomass accumulation when fine-scale variability (1-10km) is
included in the physical model. I quantified the relative contributions of biological limitations (e.g.
light limitation) and physical limitations (e.g. mixing and advection) to determine their roles in
permitting the spring bloom to occur.
4
Chapter 2
Impact of Lagrangian Sea Surface Temperature Variability on
Southern Ocean Phytoplankton Community Growth Rates
Abstract
Ocean phytoplankton play a critical role in the global carbon cycle, contributing ≈50% of global
photosynthesis. As planktonic organisms, phytoplankton encounter significant environmental variability as they are advected throughout the ocean. How this variability impacts phytoplankton
growth rates and population dynamics remains unclear. Here, we systematically investigated the
impact of different rates and magnitudes of sea surface temperature (SST) variability on phytoplankton community growth rates using surface drifter observations from the Southern Ocean
(>30°S) and a phenotype-based ecosystem model. Short-term SST variability (<7 days) had a
minimal impact on phytoplankton community growth rates. Moderate SST changes of 3-5°C over
7-21 days produced a large time lag between the temperature change and the biological response.
The impact of SST variability on community growth rates was nonlinear and a function of the rate
and magnitude of change. Additionally, the nature of variability generated in a Lagrangian reference frame (following trajectories of surface water parcels) was larger than that within an Eulerian
reference frame (fixed point), which initiated different phytoplankton responses between the two
reference frames. Finally, we found that these dynamics were not captured by the Eppley growth
model commonly used in global biogeochemical models and resulted in an overestimation of community growth rates, particularly in dynamic, strong frontal regions of the Southern Ocean. This
work demonstrates that the timescale for environmental selection (community replacement) is a
critical factor in determining community composition and takes a first step towards including the
impact of variability and biological response times into biogeochemical models.
5
Plain Language Summary
Ocean phytoplankton are fundamental to the global carbon cycle. However, it remains unclear how
environmental variability impacts phytoplankton growth, and thus, the global carbon cycle. Phytoplankton encounter environmental variability (e.g., sea surface temperature (SST) changes) as they
are transported throughout the oceans by surface currents. Here, we quantified this variability (i.e.,
in a Lagrangian reference frame) using surface drifters and investigated the impact of this variability on phytoplankton community growth rates using an ecosystem model. We also compared the
Lagrangian SST to the SST variability of a fixed point (e.g., a buoy) where ocean currents flow
past (i.e., the Eulerian reference frame) using high-resolution satellite data. We found larger SST
changes in the Lagrangian than in the Eulerian reference frame and discovered that this difference
impacted phytoplankton community structure and growth rates. The impact of SST variability was
not captured by the growth model that is typically used by global biogeochemical models. Our results provide an important extension on the classic principle that “everything is everywhere: but the
environment selects” (Hutchinson, 1961). Even when ‘everything is everywhere’, we show that the
timescale for environmental selection (community replacement) is a critical factor in determining
community composition
2.1 Introduction
Phytoplankton are a fundamental part of the global carbon cycle accounting for nearly 50% of all
photosynthesis globally (Falkowski et al., 2008). Phytoplankton also serve as the base of the marine food web and drive the ocean biological carbon pump, which acts to sequester carbon in the
deep ocean. Understanding the impact of rising global temperatures on phytoplankton communities is therefore critical for predicting the influence of anthropogenic warming on ocean ecosystems
and the global carbon cycle (Doney, 1999; Quer´ e et al., 2005). Currently, the parameterization of ´
temperature-dependent growth rates is one of the main sources of uncertainty for future carbon
cycle predictions among global biogeochemical models (Laufkotter et al., 2015). Due to their ¨
6
planktonic nature, phytoplankton will encounter anthropogenic warming in two ways: 1) as a general warming overlain on top of significant temperature variability due to advection; and 2) as
changes in variability driven by large-scale shifts in ocean physics (Boyd et al., 2016; Fu et al.,
2016; Lomas et al., 2010). Both of these processes will shift the type and magnitude of temperature
variability experienced by phytoplankton. Therefore, an improved understanding of the impact of
temperature variability on phytoplankton growth rates is necessary in order to mechanistically incorporate phytoplankton growth dynamics into biogeochemical ecosystem models and to generate
robust predictions of future changes.
Accurately assessing the type of temperature variability (rate and magnitude of change) encountered by phytoplankton in the ocean requires the correct reference frame. For phytoplankton,
the correct reference frame is Lagrangian (along trajectory) rather than an Eulerian (fixed location) reference frame (Figure 1). Modeling studies have demonstrated that both mean conditions
and variability (magnitude and rate of change) can differ markedly between the two reference
frames (e.g. Doblin et al., 2016). Here we quantify the nature of the variability experienced by
phytoplankton along Lagrangian trajectories using in situ data from the Southern Ocean – a region where global climate models lack a consensus on the impact of anthropogenic warming on
phytoplankton growth (Bopp et al., 2013). This analysis allows us to constrain realistic rates and
magnitudes of temperature changes experienced by Southern Ocean communities and determine
the impact of this variability on population growth rates.
The response of a phytoplankton community to changes in temperature is driven by individual
phytoplankter dynamics. Growth rate as a function of temperature (reaction norm) for an individual phytoplankter is unimodal and tends to be asymmetric, often with skewed tails towards
lower temperatures (Boyd, 2019). Therefore, the growth response for an individual phytoplankter
to a change in sea surface temperature (SST) depends on the starting SST relative to the optimum growth temperature (Topt, the temperature with the highest growth rate) and whether the
SST change is increasing or decreasing (Figure 2). The rate of change in growth rate will depend
on the acclimation rate (how fast the phytoplankter adjusts to the new temperature) and type of
7
acclimation of the phytoplankter (Kremer et al., 2018). When SST changes are slower than the
phytoplankter acclimation rate, the instantaneous growth rate will be equivalent to the acclimated
growth rate (i.e., the phytoplankter is able to keep up with the rate of temperature change). When
the rate of SST change is faster than the rate of acclimation, the instantaneous growth rates could
be higher or lower than the acclimated growth rate, depending on the type of response, detrimental
or beneficial, respectively (Kremer et al., 2018).
Laboratory based experiments on the impact of temperature variability on phytoplankton growth
have produced conflicting results. Some studies found an overall decrease in growth rates in a thermally variable environment relative to a stable environment (Bernhardt et al., 2018; Qu et al., 2019;
Wang et al., 2019), while others found higher growth rates under variable conditions (Schaum et
al., 2018), and some found that thermal variability did not impact community growth rates (Kling
et al., 2019; Qu et al., 2019). The lack of consensus concerning the impact of variability on phytoplankton growth rates may be due to the different magnitudes and rates of change used by the
different studies, which ranged from ≈1.5 °C/day (Schaum et al., 2018) to as high as 10 °C/day
(Bernhardt et al., 2018).
Understanding how an in situ population of phytoplankton will respond to temperature fluctuations is further complicated by phenotype and strain diversity. Multiple phenotypes can co-occur
within a population of phytoplankton each with different optimal temperatures (Topt, e.g., Webb
et al., 2009). As such, the temperature response of a population is often modeled using an Eppley curve ( αe
bT , (Eppley, 1972)) where growth rate increases exponentially with temperature
rather than as a unimodal relationship (Bopp et al., 2013). In essence, this assumes rapid phenotypic shifts within the community such that, as the temperature changes, the community rapidly
shifts its optimal growth temperature (Fig. 3b). Previous work has demonstrated that representing phytoplankton growth using an Eppley curve results in an over-estimation of phytoplankton
community growth rates (Moisan et al., 2002). In addition, the advection of communities across
large temperature gradients, such as those along a western boundary current, can result in considerable differences between the optimum growth temperature (Topt) for the community and the in situ
8
temperature (Hellweger et al., 2016). Here we build upon these studies and use a model to assess
phenotypic shifts within a population in response to different types of temperature fluctuations and
the resulting impact on population (or community) level growth rates.
In this study, we systematically assessed the effect of different magnitudes and rates of change
of temperature on phytoplankton community growth rates in the Southern Ocean (south of 30°S)
using in situ SST data and a numerical ecosystem model. This southern hemisphere region encompasses some of the lowest (0.2°C) and highest (1.6-2.0°C ) long-term mean SST variability globally
(Deser et al., 2010; Maheshwari et al., 2013). We found that relatively small changes (< 2 °C over
7-90 days) did not substantially impact community growth rates and that moderate changes (4-6°C
over 21-45 days) had the largest and longest lasting effect on community growth rates. These
moderate changes resulted in a temporary decrease in community growth rate, that lasted up to
20 generations, as the community responded to the new temperature. The response of community
growth rate to variable temperatures was non-linear and so could not easily be accounted for with
an adjustment to the Eppley curve. Finally, we found that the impact of temperature variability
on phytoplankton community growth rates was present everywhere in Southern Ocean with the
largest impact occurring in regions dominated by meso- and sub-mesoscale activity.
9
Latitude
Latitude
Latitude
Longitude Longitude
Longitude
t
1
t
2
t
3
t
1 t
2 t
3
Temperature
Lagrangian trajectories
Eulerian fixed-point
a b
c d
Figure 2.1: Lagrangian versus Eulerian reference frames. Lagrangrian reference frames follow
the water parcel itself through time. The Eulerian reference frame refers to a fixed point in space
(e.g. buoys or mooring stations) where advection of water parcels floating past the fixed point
generates temporal variability. Panels (a), (b), and (c) depict three different water masses (grey,
blue, and green) as they each pass through the fixed Eulerian location (red dot) at times t1, t2,
and t3¬. Panel (d) shows the temperature of each water mass through time (grey, blue, and green
lines) as well as the temperature recorded at the Eulerian location (red line). Note in panel (d) that
overall temperature variability in the Lagrangian reference frame (grey, blue and green lines) is
much greater than that in the Eulerian reference frame (red dots) though this may not always be
the case.
10
0 5 10 15 20 25
Temperature [o
C]
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
Growth rate [day-1]
max
Topt
Oct 27 Nov 10 Nov 24 Dec 08 Dec 22
2002
12
13
14
15
16
17
18
19
20
21
22
SST [oC]
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
0.55
0.6
0.65
Growth rate [day-1]
Topt
a
b
Figure 2.2: The impact of SST variability on individual phenotype growth rate. (a) The temperature
related growth response for a phenotype with a skewed shaped reaction norm. The values for the
optimum growth temperature (Topt) and the corresponding maximum growth rate (µmax) are shown
with dashed lines. (b) The 90-day SST profile of an example drifter trajectory (black) and the
associated changes in phenotype growth rate (blue). The orange and red arrows in the top panel
indicate the change in the phenotype growth rate associated with the corresponding changes in
SST in the bottom panel.
11
Figure 2.3: a) Map of all 90-day drifter trajectories (n = 2190) colored by SST. Two example
trajectories are highlighted in purple and magenta. b) Reaction norms for each of the 319 phenotypes in the ecosystem model. The grey lines represent all the phenotype reaction norms and the
green lines are example phenotypes to highlight the reaction norm shape. c and d) Example trajectories and their resulting model outputs. The top panels show the SST (colors), the community
growth rate estimated using the Eppley curve (dashed line), and the community growth rate from
our phenotype-based model as calculated using Eq. 5 (solid line). The bottom panel shows the
biomass through time of each phenotype (grey lines). The blue line follows the phenotype with the
highest initial biomass, the red dashed line follows the phenotype that has the highest biomass at
the end of the 90 days, and the green line follows the phenotype that has a Topt equal to the mean
SST of the trajectory.
12
2.2 Methods
The impact of SST variability on phytoplankton community growth rates was studied by combining SST observations, both in situ and from remote sensing products, and a phenotype-based
ecosystem model. Here, we focused on the impact of mixed layer SST variability on phytoplankton community growth rates and therefore did not consider growth limitations due to other sources
of variability such as nutrients, light, and mixed layer depth (e.g., (Rohr et al., 2020a; Rohr et al.,
2020b)). We tested the impact of co-limitation by temperature and nitrate and found that the results were generally consistent with the findings presented here (Appendix A1.6). Further work is
needed to investigate the impact of multiple uncorrelated environmental drivers.
2.2.1 Southern Ocean drifter profiles
Lagrangian SST data were obtained from 422 Southern Ocean surface drifters from the Global
Drifter Program with 6-hourly SST data. Float data south of 30°S from July 1999 – April 2016
was downloaded from the Drifter Data Centre at the Atlantic Oceanographic and Meteorological
Laboratory (accessed 11/2018). The lifetime of the drifters ranged from 91 days to 5.8 years with
a median duration of 521 days. Each drifter was segmented into 90-day trajectories to provide
consistency in the dataset. We used only segments that had less than 10% of missing data. This
resulted in 2,190 90-day trajectories (Figure 3a).
To estimate the magnitude of Lagrangian variability in our study region, we calculated the
range of SSTs (∆SSTmax) and the time (∆tmax) over which the temperature change occurred using
moving windows of 1 to 90 days (in 1-day increments). We then assessed the distribution of
variability across different window sizes by aggregating the data into 1°C bins for ∆SSTmax % and
1 day ∆tmax bins. For example, a 2.4°C change that occurred over 14 days was recorded in the
2-3°C and 14-day bin. To investigate the potential impact of small-scale noise, we also created
smoothed splines of each of the 90-day SST profiles using a cubic smoothing spline (csaps in
Matlab with a smoothing parameter of 0.00001). The splines filter out 25% of the variability on
13
a 1-day timescale up to 95% at the 90-day window (Figure A1.1). We then repeated the ∆SSTmax
and ∆tmax analysis on the spline data.
2.2.2 Remote sensing SST
To compare the SST variability in the Lagrangian reference frame to the variability that would be
captured in the Eulerian reference frame, we used high-resolution (0.01°horizontal resolution and
1-day temporal resolution) satellite SST data from GHRSST Level 4 MUR Global Foundation Sea
Surface Temperature Analysis (v4.1) (JPL MUR MEaSUREs Project, 2015; accessed Oct. 2018).
This dataset spanned 2003-2014 which overlaps with 71% of our 90-day drifter segments. For
each 90-day drifter segment between 2003-2014, we extracted 90 days of satellite SST data for the
latitude and longitude of the final location of the drifter, where the 90 days corresponded to the
dates of the drifter segment. We then also extracted the satellite SST along the drifter trajectories
to provide a direct comparison between the Eulerian and Lagrangian reference frames in terms of
the temporal and spatial resolution of the datasets. We performed the same ∆SSTmax and ∆tmax
variability analyses for the satellite data as the surface drifter trajectories (described in Section
2.1).
2.2.3 Idealized SST trajectories
We complemented the observed SST trajectories with idealized SST trajectories to mechanistically understand the impact of the rate and magnitude of SST change on community growth rates.
Specifically, a suite of trajectories (N = 64) was generated with both increasing and decreasing
SST trends ranging from ∆ SST = 2°to ∆ SST = 9 °C (in increments of 1 °C) over 7, 21, 45,
and 90 days. These ∆ SST values and durations were chosen based on our Lagrangian variability
analysis. To minimize initialization bias, SST was held constant for the first 30 days before increasing/decreasing. After the SST change, the SST was again held constant until the 200th day.
The final temperature for all idealized trajectories was 15 °C. The impact of the final temperature
on the model results was analyzed with a set of sensitivity experiments. The final SST had no
14
significant impact on the results when the results were reported in terms of the doubling time (generation), rather than absolute days as this normalized the effect of higher growth rates at warmer
temperatures (Appendix A1.1). Generation time was calculated as ln(2/µSS), where µSS is the stabilized community growth rate. With a final SST of 15 °C, the generation time was approximately
1.37 days.
2.2.4 Phenotype-based Ecosystem Model
To estimate the impact of variable temperature on phytoplankton community growth rates, we used
a phenotype-based ecosystem model. The model consisted of 319 phytoplankton phenotypes that
were identical in all aspects (i.e. model parameters) other than the optimal growth temperature
(Topt). Temperature dependent growth rate (µ, day−1
) was defined as a function of T (°C) (Thomas
et al., 2012):
µ(T) = α e
bT
1−
2(T −Topt)
w
2
(2.1)
where Topt was the optimal growth temperature. The value of b controlled the shape of the
reaction norm, a(day−1
)scaled the reaction norm, and w (°C) defined the width of the reaction
norm (the difference between the maximum (Tmax) and minimum (Tmin) growth temperatures). We
ran the model with two sets of reaction norms: a symmetrical, or broad, curve where b = 0 (°C-1)
and a skewed reaction norm where b = 0.3 (°C-1). Both reaction norms had a width of 14 °C (w =
20 °C), consistent with observed reaction norms for many polar species (Boyd, 2019). Sensitivity
tests were performed with reaction norm widths of 10.5 °C (w = 15 °C) and 20.5 °C (w = 29
°C) (Appendix A1.2). The results from these sensitivity tests did not differ substantially from the
simulations with a reaction norm width of 14 °C.
The parameter a scaled the reaction norms at Topt to the Eppley curve (Eppley, 1972) where
maximum growth rates ranged between 0.28 day−1
at -1.8 °C to 1.0 day−1
at 30 °C, consistent
with experimental data (Boyd, 2019). Specifically, ai was defined for each phenotype i as:
15
ai = 0.02963e
0.0405Topt (2.2)
This resulted in an increase of ≈1.5x in growth rate for every 10 degrees (i.e. a Q10 relationship
of 1.5, see Discussion). We generated 319 phenotype curves for both the broad and skewed reaction
norms with Topt ranging from -1.8 °C to 30 °C increasing by 0.1 °C (Figure 3a).
The biomass of each phytoplankton phenotype Pi was calculated at each time-step as the integral of:
dPi
dt
= µ(T)Pi −m(T)P
2
i
(2.3)
where µi(T) was the temperature-dependent growth rate for phenotype i from Equation 1. m(T)
was the temperature-dependent quadratic mortality rate
f racm3mmolCday
where:
m(T) = 0.35 ∗ a (2.4)
Here we used the same temperature dependent Eppley curve (Eq. 2) to scale mortality with temperature using SST instead of Topt where a = 1 day−1
for SST = 30 °C . We imposed a minimum
biomass (0.001 mmol C m
−3
) so that no phenotype went locally extinct, akin to the “everything
is everywhere” principle (Hutchinson, 1961). Sensitivity tests were performed with the minimum
biomass set to 0.0001 mmol C m
−3
. The minimum biomass threshold did not affect the overall
patterns but did increase both the magnitude of the difference from the community growth rates
obtained using the Eppley growth model and the time to acclimation (memory length, Section 3.2)
for both broad and skewed reaction norms (Appendix A1.3). Imposing this minimum biomass
purposefully introduced mass into the system which was accounted for by adjusting the biomass
of each phenotype to keep the total community biomass at the concentration it would have been
without the minimum biomass criteria. Specifically, the total change in biomass without the minimum biomass phenotypes was calculated using the biomass weighted community growth rate (λ)
in place of µ(T) in Eq. 3, where λ was defined as:
1
λ = ∑µi,t
Pi,t
P
(2.5)
where, µi,t was the growth rate of the ith phenotype at time t for all phenotypes with biomass
greater than the minimum, Pi,t was the biomass of the ith phenotype whose biomass was greater
than the minimum at time t, P was the sum of the biomass of all phenotypes with biomass greater
than the minimum at time t.
Several different models for mortality and grazing were tested including linear mortality, constant mortality, a dynamic zooplankton population, and a simple ecosystem model with constant
grazing pressure (see section A1.4 in Appendix). All model versions resulted in qualitatively similar results which demonstrated that the community dynamics were not particularly sensitive to
the top-down control formulation in the model (Section A1.4 in Appendix). Here, we present the
quadratic mortality as it was the simplest model with smooth (non-oscillatory) solutions.
The ecosystem model was forced with each of the 2,190 drifter segments (see Figure 3c-e for
examples), the corresponding smoothed splines, the idealized SST trajectories, and the satellitederived SSTs. The initial biomass of phenotypes with a Topt within ±2.5 °C of the starting SST
value were randomized to simulate previously accumulated biomass with phenotypes outside this
range set to the minimum biomass. Simulations that used idealized SST trajectories were performed 100 times with different initial biomass conditions to account for stochasticity in the model
initialization.
2.2.5 Acclimation Rate
To test the impact of different acclimation timescales, we performed sensitivity tests in which we
incorporated a linear acclimation rate for all phenotypes in the model. Specifically, we incorporated
a timescale over which an individual phenotype could change its growth rate in response to a
temperature change. For example, if SST rapidly changed from 15 °C to 16 °C, a phenotype with
an acclimation timescale of 0.2 °C day−1 would move from the growth rate at 15 °C to the growth
rate at 15.2 °C in one day. If the SST then held constant at 16 °C, the phenotype would acclimate
17
by the end of the fifth day. We tested acclimation rates ranging from 0.2 °C day−1
to 0.6 °C
day−1
in increments of 0.1 °C day−1 which are consistent with acclimation rates determined for
the Southern Ocean diatom F. cylindrus (see Table A1.2).The model with acclimation was forced
with the idealized SST trajectories for a ∆SST = 2 °C in 7 days (0.29 °C day−1
), 3 °C in 7 days
(0.43 °C day−1
), 4 °C in 7 days (0.57 °C day−1
) and 5 °C in 21 days (0.24 °C day−1
). These
intervals corresponded to the magnitudes and rates of change most commonly experienced by the
drifter trajectories (see Section 3.1) for which the rate of change was greater than 0.2 °C day−1
.
2.3 Results
2.3.1 SST variability
We characterized in situ SST variability experienced by phytoplankton (i.e., in a Lagrangian reference frame) using the surface drifter SST data. Seasonal dynamics were not filtered out as
they were important sources of SST variability encountered by phytoplankton. While the surface
drifters may have been subjected to some physical movements that phytoplankton do not encounter
(e.g., lateral transfer across fronts due to wind rather than subduction and mixing), they provided
the best in situ dataset for studying Lagrangian variability in surface temperature. However, to
minimize the impact of unrealistic fluctuations in the drifter dataset, we limited our subsequent
analyses to the most frequently measured scales of variability within the drifter record. The average ∆SSTmax values ranged from 0.9 °C ±0.7 °C (1σ) for the 7-day window, which corresponded
to 0.13 °C/day change over the 7 days, to 4.2 °C ±2.0 °C (1σ) for the 90-day window or 0.05
°C/day change (Figure A1.12, Table A1.1). The latter was consistent with the expected seasonal
SST cycle for the Southern Ocean (Reynolds et al., 1994). The SST variability of the drifters (standard deviation over the window) was highly correlated with ∆SSTmax (R2 = 0.92, p ¡ 0.01, Figure
A1.13).
Using the ∆SSTmax analysis, we were able to quantify the most common types of variability
encountered in situ in terms of both the magnitude of change and the rate of change (Figure 4).
18
Due to the difference in the number of data points generated by the moving windows, we assessed
the frequency of each ∆SSTmax within a given window length (y-axis) such that the highest value
across the row indicates the most likely ∆SSTmax for that window length. The ∆SSTmax bins sum
to 100% across the row. Overall, there is a trend of increasing ∆SSTmax with increasing window
length, as expected. We selected four representative window lengths, 7 days, 21 days, 45 days, and
90 days, as illustrative examples though the results are not dependent on these selections. A 7-day
window was most likely to have a ∆SSTmax of 2 °C or less (82%), and ≈3% of the trajectories
recorded a ∆SSTmax of 4 °C. Over a 21-day window, most trajectories had a ∆SSTmax of 2-3 °C
(combined accounting for 86% of data) and ≈10% of the trajectories had a ∆SSTmax of 4-5 °C.
∆SSTmax reached as high as 9 °C for the 90-day windows but accounted for only 2.5% of the data
in that window.
A comparison of Lagrangian and Eulerian reference frames demonstrated that, while the overall
patterns of variability were similar, the Lagrangian reference frame was more likely to capture
large ∆SSTmax (Figure A1.14). This was true for both the drifter and satellite derived Lagrangian
trajectories when compared to the satellite SST in the Eulerian reference frame. For example,
within a time-frame of 21 to 30 days, a ∆SSTmax greater than 3 °C was more likely to occur in both
the drifter and satellite derived Lagrangian (17%) trajectories than in the Eulerian (11%). Similarly,
for the 51-day to 60-day windows, ∆SSTmax of 2-4 °C were common in both reference frames, but
changes ¿4 °C were more common in the satellite derived Lagrangian trajectories (24%) and the
drifter Lagrangrian trajectories (23%) than the satellite Eulerian data (16%). This same pattern
was consistently observed for all windows from 1 to 90 days. The impact of these differences in
SST changes on phytoplankton community growth rates are discussed below (see Sections 3.2 and
3.3).
For most of the SST data recorded by the drifters, the rate of SST change was slower than
the expected phytoplankton acclimation rates. Acclimation rates for the Southern Ocean diatom
F. cylindrus are on average 0.3 °C/day (Table A1.2). For the drifter trajectories, only 8% of the
1-day bins (n = 197,100 days) recorded rates of SST change greater than 0.3 °C/day and less than
19
2% of the daily bins recorded rates of change greater than 0.6 °C/day (Figure A1.15).Because SST
rates of change were typically slower than the phytoplankton acclimation rate, we hypothesize
that, for the majority of the Southern Ocean, the rate of acclimation will not play a major role in
the community response. Therefore, to simplify model dynamics, we ran our model with rapid
acclimation such that each phenotype responded directly to SST changes. See Section 4.1 for
discussion about situations in which acclimation may be important.
2.3.2 Impact of variable SSTs on community growth rates
We used idealized simulations to develop a mechanistic understanding of how variability impacts
community growth rates. For small, gradual SST changes of less than 2-3 °C in 45-90 days (0.02-
0.07 °C/day), the community growth rates changed linearly with the SST changes during the period
of SST transition and then stabilized once SST stopped changing. When the rate of change was
slow, the distribution of phenotypes within the community changed at the same rate as the SST such
that the Topt of the most dominate phenotype closely matched the SST. As a result, the temporal
response in the community growth rate from the phenotype model was similar to the growth rate
from a null model using an Eppley curve parameterization.
For SST rates of change larger than 2-3 °C in 45-90 days (0.02-0.07 °C/day), community
growth rates initially increased or decreased depending on the sign of the SST change, but then
began to decrease rapidly (see Figure A1.16 for example). Once SSTs stabilized at the final value,
community growth rates increased and eventually stabilized. Out of the environmentally relevant
rates of SST change, 4 °C in 7 days (0.57 °C/day) resulted in the largest change 70% ±1% (1σ)
during the low growth period (Figure 5a). While the absolute percent change in growth rate was
sensitive to model formulation and parameter values, the qualitative relationships presented here
were robust (Figure A1.8).
The impact of temperature variability on community growth rates is a function of both changes
in the growth rates of individual phenotypes (i.e. shifts along a reaction norm) and shifts in the
community composition (i.e. abundance of different phenotypes). The low growth phase after
20
a shift in SST (either increasing or decreasing) was caused by the SSTs extending beyond the
thermal optimum of the initial community such that the bulk of the biomass was growing slowly.
During this period, the individual phenotypes with elevated growth rates only made up a small
fraction of the community and so did not contribute significantly to the community growth rate.
The community growth rates then rebounded as these high growth phenotypes increased their
biomass and eventually became the dominate biomass group. Faster rates of SST change moved
the community out of the thermal optimum of the initial community more quickly than smaller
rates of change, and therefore larger and faster ∆SSTs resulted in greater decreases in community
growth rates. However, the high growth individuals were able to dominate the community more
quickly due to the high loss rates for the slow (or no) growth individuals and so the community
growth rates rebounded more quickly for rapid relative to moderate rates of SST change. For rapid
SST changes, the rate and type of acclimation response could potentially play a role in the shifts
in community growth rates depending on the nature of the plastic response (see Discussion).
An Eppley curve was unable to capture the impact of variability in SST on community growth
rates due to the non-linear phenotype dynamics. Community growth rates derived directly using
the Eppley curve model were always larger than those simulated by the phenotype model, consistent with previous work (Moisan et al., 2002; Bernhardt et al., 2018). The difference between the
phenotype modeled growth rates and the Eppley curve estimates varied as a function of SST variability (Figures 5b, S18). As the ∆SST increased over a given window length, so did the difference
between the phenotype model and the Eppley curve estimate. The largest departures occurred for
∆SSTs of 4 °C and 5 °C over 7 and 21 days, respectively, with up to 80% lower simulated community growth rates for the phenotype model. Generally, larger ∆SSTs and faster rates of change
(changes occurring over a few generations) resulted in larger differences between the models.
Although we focus on the impact of temperature-limitation on phytoplankton growth in this
study, nutrient limitation also plays an important role in co-limiting phytoplankton growth in the
Southern Ocean (Cochlan, 2008). While a full analysis of the impact of fluctuating co-limitation
is beyond the scope of this study, we conducted a set of model simulations to test the impact of
21
temperature and nitrate limitation on the observed dynamics (Appendix A1.7). Given the nonlinear formulation of nutrient limitation and the relatively low half-saturation values for the uptake
of nitrate, for the majority of the Southern Ocean the model results with co-limitation are similar to
those from simulations with temperature limitation only. This is because the variation in nutrient
limitation was small compared to the variation in temperature limitation (at the level of individual
phenotypes).
2.3.3 Memory Effect
The timescales of the biological response to temperature fluctuations varied as a function of the
overall magnitude and direction (increasing or decreasing) of SST change, the duration of the SST
change, and shape of the reaction norm (broad vs skewed) for the individuals within the population.
Here we define the timescale of biological response as the “memory effect” – the time for the
community growth rate to stabilize (±5% of the stable value). Here we present the memory effect
in terms of generations calculated using the final stable growth rate. This allowed us to understand
the relative impact of temperature change on phytoplankton using a common currency such that
our results are not growth rate dependent.
The most common ∆SST changes (Figure 4) were associated with the longest memory effects
(Figures 5b, S19). Nearly all of the environmentally relevant ∆SST values were sufficient to create
a memory effect of longer than 2 generations. Moderate changes of 3-4 °C over 7-45 days or 4-28
generations (0.07-0.57 °C/day) resulted in the longest memory effects of up to up to 22 generations
for both reaction norm shapes (Figure 5c, Figure A1.19). This biological response time is nearly
five times longer than the duration of the temperature change. Larger SST changes (5-6 °C) that
occurred over 45 or 90 days or 28-57 generations (0.06-0.13 °C/day) tended to have shorter memory effects (≈8-19 generations) than moderate changes that occurred over the same time frame, but
this difference was not statistically significant. Longer memory effects for moderate SST changes
resulted from dominant phenotypes in the previously acclimated community being able to grow in
the new environment, albeit at a reduced rate. This increased the time required for the phenotypes
22
optimally suited for the new environment to dominate the community, which resulted in larger
memory effects.
The sign of the SST change also impacted the response time of the community. Decreasing
∆SSTs had longer memory effects by an average of 7 generations compared to increasing ∆SSTs
(t-test, 95% CI) for the skewed shaped reaction norms. The longer memory effect was due to the
long tail on the decreasing side of the reaction norm, which allowed the phenotypes in the initial
community to grow during decreasing SST conditions (Figure A1.19). For reaction norms that
were symmetrical about the optimum growth temperature, the direction of ∆SST did not matter,
and the memory lengths were not statistically different for increasing and decreasing ∆SSTs (t-test,
95% CI) (Figure A1.19).
2.3.4 Southern Ocean Drifter Trajectories
The idealized simulations allowed for a mechanistic characterization of how phytoplankton community growth rates vary as a function of rate and magnitude of SST change (sections 3.2 and 3.3).
However, in the ocean, SST change is much more complicated as phytoplankton are exposed to a
large variety of rates and durations of SST changes. We used Southern Ocean drifter trajectories
to investigate the impact of in situ SST variability on community growth rates. When phenotypic
diversity was considered (phenotype model), variable SST resulted in lower average community
growth rates compared to the Eppley curve approximation (Figure 6). Though nutrient and light
limitation were not explicitly included in the phenotype model, the simulated growth rates were
consistent with in situ (Buitenhuis et al., 2013), remote-sensed based (Arteaga et al., 2020), and
incubation derived growth rates (Boyd, 2019; Boyd et al., 2013), and therefore reasonably captured
growth dynamics. As drifter SST variability increased, so did the difference between the phenotype model and the Eppley curve approximation, consistent with the idealized simulation results.
The mean percent difference between the phenotype model and the Eppley curve approximation
over the 90-day trajectories ranged between -142% to -11.5% with a mean of -25.8% (±16.6%
1σ) for the skewed reaction norms. A similar pattern was observed for the broad shaped reaction
23
norms, but the magnitude of the difference was smaller and ranged from just -1.3% to -53.2% different with a mean of -6.1% (±5% 1σ) (Figure A1.20). Trajectories with higher mean SSTs were
affected less by SST variability than trajectories with lower SSTs because faster growth rates at
higher temperatures allowed quicker responses to SST changes.
To isolate the impact that short-term variability may have on community growth rates relative
to longer-term shifts, we compared the 90-day mean biomass-weighted community growth rate of
the drifter trajectories to the smoothed splines derived from the trajectories. Removing short-term
variability had no significant impact on community growth rates (t-test, 95% CI; Figure A1.21).
24
10
20
30
40
50
60
70
80
90
Window Length [days]
Drifter ∆SSTmax
0
5
10
15
20
25
% of Data
∆SSTmax [o
C]
0 1 2 3 4 5 6 7 8 9
Figure 2.4: SST variability analysis. The frequency of ∆SSTmax changes from the drifter segments
over different window lengths are shown. Data are presented as total percent of data that fall within
that window length such that each row sums to 100%. There is a general pattern of increasing
magnitudes of ∆SSTmax over longer window lengths.
25
∆SST
23456789
-60
-50
-40
-30
-20
-10
0
7
21
45
90
% Decrease
23456789
7
21
45
90
0
5
10
15
20
Change in community growth rate
23456789
∆SST
7
21
45
90
-70
-60
-50
-40
-30
-20
-10
0
% Difference [generations]
∆SST
Difference from Eppley Curve
Memory Length
a
b
c
Figure 2.5: Simulated response of phytoplankton community with skewed shaped reaction norm to
increasing ∆SST (see Appendix for decreasing ∆SST conditions and broad reaction norm results).
Panel (a) plots the decline in community growth rate in the phenotype model that results from the
SST moving out of the thermal niche of the original population (see Methods and Figure A1.16).
Data that are greyed out represent ∆SST and window length combinations that were not supported
by the results from Figure 3. Panel (b) shows the percent difference between the Eppley growth
model approximation and the phenotype modeled community growth rates at the point where SST
stabilizes (see Figure A1.16 for example). Panel (c) plots the memory effect length associated
with SST changes in the idealized simulations. This represents the time it takes for the community
growth rate to be within 5% of the steady state growth rate at the final SST from the first time-step
that SST is constant (See Figure A1.16 for example).
26
SST standard deviation [1-σ]
-60
-55
-50
-45
-40
-35
-30
-25
-20
-15
-10
Mean % difference over 90-day trajectory
Percent reduction in mean community growth rate
relative to Eppley curve
0
5
10
15
20
25
Mean SST [oC]
Trajectory
Idealized
Example
0 0.5 1 1.5 2 2.5 3
Figure 2.6: Impact of SST variability on community growth rate. The average percent difference
in community growth rate between the phenotype model and the Eppley growth model from the
90-day drifter segments are plotted against the standard deviation (1σ) of the drifter SST. Each
segment is colored by the mean SST. Results from the idealized trajectories are shown as black
diamonds with filled circles denoting increasing SST trajectories and open circles denoting decreasing SSTs. Pink triangles represent the two example trajectories from Figure 3. Results shown
here are for skewed shaped reaction norms, see Figure A1.20 for results for the broad shaped reaction norms.
27
2.4 Discussion
2.4.1 Impact of Acclimation
As ocean surface temperature shifts, two processes occur simultaneously: 1) individual phytoplankton phenotypes respond to the change in temperature (acclimation), and 2) phenotype abundance within the community shifts towards individuals with higher maximum growth rates at the
new temperature. In this study, we investigated the impact of these individual-level dynamics
on community-level growth. We demonstrate that shifts in phenotype abundance are the primary
drivers of community growth rate dynamics. This is in large part due to in situ rates of SST changes
being slower than the rates of individual acclimation (based on laboratory estimates), even for the
dynamic Southern Ocean. When individual acclimation rates were slower than the rate of SST
change, we observed a delay in the low growth phase and a smaller magnitude decrease in community growth rates (Figure A1.22). The memory effect increased linearly as individual acclimation
rates decreased (longer acclimation time) (Figures 7 and S23).
The representation of phenotype acclimation in the model was a simplistic representation of
phenotypic plasticity (see Methods). In reality, plastic responses are much more complex and
nonlinear and most likely vary among species (Kremer et al., 2018). Additional work is needed
to better constrain both the range of acclimation timescales and the mechanisms of phenotypic
plasticity. However, our results suggest that these dynamics will only become important under
rapid temperature changes which are infrequent in the ocean.
2.4.2 Implications for in situ community composition
Our findings support the important role of thermal history in shaping the response of phytoplankton
communities to changes in temperature. We have shown that SST variability can lower community
growth rates for tens of generations following SST perturbation. This indicates that, for many
regions of the ocean, the phytoplankton community will not be fully acclimated to local conditions
as a result of the mismatch between timescales of physical variability and biological response.
28
Acclimation Rate [o
C/day]
0
5
10
15
20
25
Memory Length [generations]
2o
C/7 days (0.29o
C/day)
3o
C/7 days (0.43o
C/day)
4o
C/7 days (0.57o
C/day)
5o
C/21 days (0.24o
C/day)
0 0.1 0.2 0.3 0.4 0.5 0.6
Figure 2.7: The impact of acclimation on the number of generations for which the effect persists
(memory length). Acclimation rates that were slower than the rate of SST change resulted in longer
memory lengths than for simulations in which acclimation rate was equal to or faster than the SST
rate of change.
This mismatch in timescales will be a function of the rate and magnitude of SST variability that
phytoplankton in the water mass were previously exposed to and may be reflected in physiological
properties such as optimum growth temperature or overall community growth rate.
Our results also provide an important extension on the classic principle that “everything is everywhere: but the environment selects” Hutchinson, 1961. Even when ‘everything is everywhere’,
we show that the timescale for environmental selection (community replacement) is a critical factor
in determining community composition. Specifically, we hypothesize that even when the ‘optimal’
organism is present in an environment, environmental variability generated by local physics, lateral
advection, and seasonal trends can delay or prevent that organism from dominating the community.
29
This hypothesis is supported by previous modeling work that has shown a time-lag on the order
of weeks to a month in the phytoplankton community growth response to SST changes due to
lateral advection and seasonal trends (Moisan et al., 2002; Hellweger et al., 2016). Here, we have
quantified the relationship between varying rates of SST variability and the timescale required for
community replacement to impact the community composition.
We tested the impact of Eulerian versus Lagrangian variability on community growth rates
and demonstrated significant differences for locations in which SST variability differed in the two
reference frames. Specifically, while the final SST of the drifter segments and satellite data were
not statistically different (t-test, 95% CI, Figure A1.24), differences in the nature of variability in
the proceeding 90 days resulted in a significant difference between the final SST and the Topt of the
most abundant phenotype (t-test, 95% CI, Figures 8 and S25). The magnitude of the offset between
SST and Topt depended on the timing of SST changes throughout the 90-day profiles. When SST
changes were slow, the offset between SST and the Topt of the most abundant phenotype were
negligible (Figure A1.26 for an example). Large SST changes that occurred early in the 90-day
segment allowed sufficient time for the community to respond (e.g. Figure A1.26). When SST
changes occurred later in the 90-days, the community did not have sufficient time to respond
which caused a larger offset between the SST at day 90 and the Topt of the community (e.g., Figure
A1.27). Different phenotype distributions for the Eulerian versus Lagrangian reference frames is
consistent with previous results that showed advection of phytoplankton communities was a key
process in shaping phytoplankton diversity (Barton et al., 2010; Clayton et al., 2013; Levy et al., ´
2014).
The shape of the reaction norm impacts the community response to temperature variability
and phenotype competitive advantage. Under decreasing temperatures, a phenotype with a skewed
reaction norm (Topt closer to Tmax than Tmin) has a competitive advantage over a phenotype with
a broad reaction norm (Topt at the center of Tmax and Tmin), given the same reaction norm width
and Topt. A skewed reaction norm provides a larger range of temperatures < Topt under which
the phenotype can grow. Therefore, organisms with skewed reaction norms should be adapted to
30
have Topt values close to maximum encountered temperatures not only due to the rapid decline in
growth rates for temperatures greater than Topt but also due to the competitive advantage under
temperatures less than Topt. Conversely, broad reaction norms are favored when temperatures are
warming, as expected, or when temperatures are more variable. In a highly variable region such as
the Southern Ocean, there should be selective pressure for either broad reaction norms with large
growth ranges beyond Topt (Moisan et al., 2002) or skewed reaction norms where Topt is higher
than mean SSTs (Thomas et al., 2012).
2.4.3 Implications for simulating community growth rates in global biogeochemical
models
A form of the Eppley curve Q10 temperature-growth response
µ = µoQ
T−To
10
10
is widely used in
global biogeochemical models (Bopp et al., 2013), where typical model Q10 values range between
1.5-2 (Sherman et al., 2016). The premise behind employing a Q10 growth equation is that each
modeled functional group encompasses many species or strains and so the Eppley curve may be
a reasonable representation of the group dynamics. However, as we have demonstrated, community growth rates (or functional group growth rates in the model framework) will depend on the
underlying phenotype dynamics, which are a function of the rate, magnitude, and direction of
temperature change and the shape of the species/strains’ thermal response curve. As a result, the
Q10 temperature-growth response not only underestimates temperature-limitation on community
growth rates (i.e., overestimates growth rates) but does so as a function of SST, SST variability,
and reaction norm shape. Our work indicates that adjusting the Q10 relationship to use a lower
exponent as previously suggested (Sherman et al., 2016) will only partially capture realistic dynamics. Because phytoplankton play a key role in sequestering carbon dioxide from the Earth’s
atmosphere, by overestimating phytoplankton growth rates, and thus overestimating carbon uptake,
biogeochemical models may be underestimating the extent of future anthropogenic warming.
To predict changes in phytoplankton community growth rates robustly, models must also consider the impact of different types of SST variability and the appropriate reference frame for this
31
variability. Specifically, we have shown that SST variability can differ markedly between the
Eulerian reference frame and the Lagrangian reference frame (Figure A1.14). While the spatial
patterns of SST variability in the Southern Ocean were similar between Eulerian and Lagrangian
reference frames (Figure 9a-c), the Eulerian reference frame exhibited substantially less variability. Consequentially, the offset between Eppley curve approximation and the phenotype model
was substantially less for the Eulerian relative to the Lagrangian reference frame (Figure 9d-f).
This pattern was consistent for both the drifter Lagrangrian trajectories and the satellite derived
trajectories.
Models such as DARWIN (Follows et al., 2007) resolve phenotypes with a range of thermal
reaction norms and so will capture the community growth rate dynamics presented here. However,
additional work is needed to compare the variability encountered by functional group phenotypes
in large-scale models integrated in an Eulerian framework to true Lagrangian variability.
Improving the parameterized temperature-growth relationship is particularly important in the
Southern Ocean given the uncertainty of future primary productivity in this ocean basin (Bopp et
al., 2013). We used our model results to identify key regions within the Southern Ocean that might
be most strongly impacted by temperature variability. Three particular regions stand out that exhibited the most SST variability and had the largest relative deviations from the Eppley curve (Q10)
approximation: the Malvinas-Brazil confluence zone; the Agulhas Retroflection region; and downstream from these two along the Subtropical Front near ≈45 °S, 60 °E (Figures 9a-c). All three
regions were previously identified as highly dynamic, strong frontal regions (Artana et al., 2019;
Beal et al., 2015; Graham et al., 2013) and shown to be important hot-spots for phytoplankton
diversity (Barton et al., 2010; Clayton et al., 2013; d’Ovidio et al., 2010; Soccodato et al., 2016). It
is possible that in these highly dynamic frontal regions the floats were subjected to physical movements across the fronts that was previously thought to elude phytoplankton movements. However,
recent field and modeling studies have shown that cross-front transfer and diapycnal mixing can
occur due to the fine-scale physics associated with these strong fronts (Clayton et al., 2017; Mahadevan, 2016; Wenegrat et al., 2020). Our results also showed that large SST changes were not
32
required for temperature variations to have a lasting impact on community growth rates. Regions
of the Southern Ocean with moderate (1-2 °C, 1σ) SST variability also recorded equally large differences in community growth rate, often at least 30% smaller than Eppley curve approximations
and up to >=80% smaller than Eppley curve approximations.
33
-6 -4 -2 0246
Drifter Offset [o
C]
-6
-4
-2
0
2
4
6
Satellite Offset [oC]
Figure 2.8: The impact of Lagrangian and Eulerian variability on community composition. Here
we plot the difference between the Topt of the most abundant phenotype at the end of each 90-day
trajectory and the final SST for the drifter trajectory (x-axis) and the satellite data (y-axis). The
final SSTs for the drifter and satellite data are not statistically different (t-test, 95% CI). Therefore,
deviations from the 1:1 line demonstrate the impact of a Lagrangian versus Eulerian reference
frame on community composition.
34
Figure 2.9: Distribution of SST variability (a-c) and the deviation in community growth rate from
the Eppley growth model (d-f) over the Southern Ocean (>30 °S). Only those drifters which overlap in space and time with the satellite data are shown. For full results, see Figure A1.29. Three key
regions of high SST variability stand out: Malvinas-Brazil confluence zone, the Agulhas Retroflection, and the Subtropical front. These regions have enhanced SST variability in all datasets but
higher variability in the Lagrangian trajectories. These high variability regions correspond to large
differences between the phenotype model growth rates and the Eppley approximation of growth, a
pattern consistent across all three sets of simulations.
35
2.5 Conclusions
We utilized idealized SST simulations and SST data from ocean surface drifters to show that synoptic SST variability on timescales of a few days to a few weeks decreases phytoplankton community growth rates, while higher frequency variability has little impact. The time taken for the
community growth rate to reflect the new environment was dependent upon the rate and magnitude of temperature change, the direction of change, and the shape of the thermal response curve.
The largest memory effects resulted from moderate changes in SST that occurred over 1-3 weeks.
This impact of SST variability can cause a large offset between a phenotype-based temperaturedependent community growth rate and an Eppley curve-based approximation and suggests that
phytoplankton communities sampled in situ may often not be adjusted to local conditions. Given
the highly variable nature of the ocean and importance of environmental variability for phytoplankton physiology, it is critical to consider the appropriate reference frame and the magnitude
and duration of variability when studying phytoplankton dynamics. Here we demonstrated that
variability captured in the Lagrangian reference frame (by drifters) was, in many instances, different from variability in the Eulerian frame and that this had significant impacts for estimating
phytoplankton growth rates. These findings have potentially far-reaching implications for how
temperature-dependent phytoplankton growth is represented in global biogeochemical models.
2.6 Acknowledgements
We would like to thank Robert Strzepek for providing us with the acclimation rate data (Table A1.2). Additionally, we acknowledge funding support from the National Science Foundation (NSF OCE 1538525 to NML). PWB was supported by an ARC Laureate fellowship and
the Australian Department of Industry funded Australian Antarctic Program Partnership. SCD
acknowledges support from NSF Office of Polar Programs (grant PLR-1440435 to the Palmer
Long Term Ecological Research project). JNH acknowledges support from the University of
Gothenburg Natural Sciences Sabbatical Program. Drifter data used here can be obtained from
36
the Drifter Data Centre at the Atlantic Oceanographic and Meteorological Laboratory (https:
//www.aoml.noaa.gov/phod/gdp/) and satellite data are available from the GHRSST Level 4
MUR Global Foundation Sea Surface Temperature Analysis (https://podaac.jpl.nasa.gov/
dataset/MUR-JPL-L4-GLOB-v4.1).
37
Chapter 3
Lagrangian mixing variablity
Abstract
Submesoscale processes play a crucial role in shaping the physical and biogeochemical dynamics
of eastern boundary upwelling systems. Accurately representing these fine-scale features in biogeochemical models is critical for understanding and predicting ecosystem responses and carbon
cycling. However, traditional Eulerian (fixed-point) measurements may not fully capture the dynamic environment experienced by organisms and particles as they are transported by currents. Using Lagrangian particle tracking in a high-resolution (500m) model of the California Current System (CCS), we quantify differences between Eulerian and Lagrangian reference frames in terms of
environmental variability (temperature, mixing, and submesoscale feature encounter). Our results
challenge the notion of consistent Eulerian overestimation of mixing, revealing nuanced spatial
relationships between Eulerian and Lagrangian variability. We observe distinct spatial patterns
in Lagrangian particle interaction with submesoscale features where nearshore particles encountered features more frequently than expected based on their occurrence. Establishing a robust
relationship between Eulerian observations and Lagrangian sampling is essential for incorporating submesoscale dynamics into biogeochemical models. This study highlights the importance of
such a connection for accurately predicting ecosystem responses and carbon cycling in the face of
a changing environment.
38
3.1 Introduction
Eastern boundary upwelling systems, like the California Current System (CCS) play a critical role
in the global carbon cycle. Seasonal upwelling brings nutrient rich water to the sun lit surface fueling primary production by phytoplankton. The conversion of carbon dioxide into organic carbon
by phytoplankton plays an important role in the ocean carbon cycle (Ricour et al., 2023). The rate
of photosynthesis and biomass accumulation, as well as phytoplankton community composition,
are in large part set by the history of environmental variability experienced by the phytoplankton
(Hellweger et al., 2016; Baudry et al., 2018; Boyd et al., 2016; Doblin et al., 2016). As phytoplankton are largely passive drifters, they encounter environmental variability in the Lagrangian
reference frame, along the water mass trajectory. To understand and correctly model phytoplankton growth, it is essential to capture the response of these communities to variability experienced
along a Lagrangian path (Zaiss et al., 2021).
Previous work has investigated Lagrangian trajectories in mesoscale resolving models (O20-
100km grid size) and demonstrated that the environmental history along a Lagrangian trajectory
differs from the Eulerian environmental history at the same location (e.g. Doblin et al., 2016).
The difference between the two reference frames has been demonstrated to cause differences in
phytoplankton growth dynamics (Zaiss et al., 2021; Baudry et al., 2018).
Physics at smaller scales that are not resolved in these models (e.g. submesoscale dynamics,
O1-10km scale) also have the potential to significantly impact phytoplankton growth by rapidly
altering the chemical and physical environment. Specifically, it remains uncharacterized how Lagrangian particles interact with features in terms of time spent in the submesoscale features and
the resulting advective transport. Therefore, it is important to quantify the fine-scale environmental
variability experienced by phytoplankton in the Lagrangian frame relative to the Eulerian.
Physical dynamics occurring at the sub-mesoscale (1-10 km) range from features such as eddies
and fronts to internal gravity waves and enhanced diffusive mixing. These features are challenging to observe in situ due to their small spatial scales and short temporal scales (O days). Both
39
high-resolution satellite observations (Morrow et al., 2019) and autonomous vehicles with highresolution sampling capabilities (Swart et al., 2020; Omand et al., 2017) are providing key datasets
but have several key limitations: 1) these datasets are confined to two dimensional dynamics (satellites), 2) they are often concentrated near convergent zones due to surface currents(surface drifters,
Ohlmann et al., 2017) or 3) they are subject to large biases depending on sampling frequency (gliders, Patmore et al., 2024). Consequently, high-resolution models are still one of the primary tools
for studying submesoscale dynamics.
Here, we investigate how physical dynamics and environmental history differs between Lagrangian and Eulerian reference frames in a high-resolution (500m) physical ocean model of the
CCS. We quantify differences in temperature and mixing between the two reference frames in addition to differences in the encounter rate of the Lagrangian particles with submesoscale eddies
and fronts relative to the abundance of these features. Characterizing the differences in these environmental variables between the two reference frames is a first step in examining how fine-scale
physics can impact the global carbon cycle.
3.2 Methods
To investigate the impact of fine-scale physical dynamics, such as submesoscale features, on Lagrangian particles, we advected particles in a high-resolution, submesoscale resolving, physical
ocean model for 66.7 days.
3.2.1 Physical Model
Here we use a high resolution, physical oceanographic, regional model simulation for the CCS.
The regional model was obtained via nesting simulation of a global ocean numerical simulation
known as LLC4320 from the Massachusetts Institute of Technology general circulation model
(Torres et al., 2018; Flexas et al., 2019). The LLC4320 simulation has a horizontal resolution of
1/48≈2km) and was run for the period from 13 September 2011 to 15 November 2012. The CCS
domain was 31Nto42N latitude and -132to − 119.5 longitude and was forced at the boundaries
40
by hourly output from the LLC4320 simulation. The CCS simulation horizontal resolution was
increased to 500m with the same 87 telescoping vertical levels as used in LLC4320. with 500m
horizontal resolution and 90 vertical levels and was forced at the boundaries by hourly output from
the LLC4320 simulation. Previous work has demonstrated that this high-resolution simulation
captures the mixing and advective properties of fine-scale dynamics such as internal gravity waves
and submesoscale eddies and fronts(Torres et al., 2022). Here we focused on a nearly 67-day
period (66.7 days) during the upwelling season from March 27 0900 to June 2 0100.
3.2.2 Lagrangian Trajectories
To understand the physical environment experience by passive drifters in the CCS, we calculated
3-dimensional Lagrangian trajectories within the model. The 3D trajectories were calculated at
hourly time steps, using linearly interpolated horizontal and vertical velocities and a first-order
Euler advection method (Van der Stocken et al., 2019). The code for calculating the Lagrangian
trajectories can be found at https://github.com/LevineLab.
A total of 2,065 Lagrangian particles were released on day 1 of the simulation at 100m depth
and the resulting trajectories were calculated over the 66.7 day simulation. Specifically, seven sets
295 Lagrangian particles spaced 1km apart from each other were released along each of seven
lines of latitude from 34 to 40 at 1 ncrements (Figure 1a). We investigated the impact of releasing
particles at different depths and found that it did not influence the primary results (Appendix A2.1).
Particles that advected close to the model borders were discarded to avoid spurious velocity values
caused by the edge of the model domain (north of 40.5 and south of 32.5 , west of -132.5 and
east of 120.5 ; Figure 1b). If a particle was trapped against coastal bathymetry or was advected
below 200m, that trajectory was also excluded. After excluding these particle trajectories, 1463
full trajectories remained and were used for subsequent analyses.
41
a b
15m
50m
100m
150m
200m
c d
42 44 46 48 50 52 54 56
Simulation Day
-80
-70
-60
-50
-40
-30
-20
-10
0
Depth [m]
-2
-1.5
-1
-0.5
0
0.5
1
1.5
2
Normalized Vorticity
-123.6 -123.2 -122.8 -122.4 35.4
35.6
35.8
36
36.2
36.4
-1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
1
Normalized vorticity
5
10
15
20
25
30
35
40
45
50
Simulation Day
SST [oC]
SST [oC]
Figure 3.1: Examples of Lagrangian particle trajectories. a) Starting locations of the Lagrangian particles at
the first timestep, colored according to starting latitude. The background is colored by sea surface temperature [] at the first time step. b) Final location of the particles on day 66. The color represents the latitude that
the particles started at as shown in (a). The size of the particle represents the depth of the particle with larger
sizes indicating deeper depths as specified in the legend (meters). c) Trajectory of an example Lagrangian
particle through time. The background is colored by the relative vorticity at day 54.5 and the trajectory is
colored by simulation day. At simulation day 54.5, the Lagrangian particle encounters a cyclonic frontal
feature, marked by a star. (d) The depth through time of the example trajectory shown in (c). The line is
colored by the relative vorticity that the Lagrangian particle encountered. The initial encounter with the
frontal feature is marked by a start at simulation day 54.5 after which the particle is rapidly subducted by
the frontal flow.
3.2.3 Identification of Submesoscale features
A primary aim of our study was to determine how often the Lagrangian particles encountered a
submesoscale eddy or front. Submesoscale features were identified based on vorticity calculated
from the model velocity fields. Specifically, we used normalized relative vorticity values (Rossby
number) which were calculated as:
ζ =
dv
dx −
du
dy
f
(3.1)
42
where du and dv are the changes in the horizontal velocities in the dy and dx directions, respectively. f is the Coriolis parameter calculated as:
f = 2Ω∗ sin(Latitude) (3.2)
where Ω is the angular velocity of the Earth and Latitude is in radians. Vertical velocities in the
CCS are dominated by internal gravity waves so a high-pass filtering method was used to remove
the movement due to internal gravity waves and isolate the velocities associated with submesoscale
features following the procedure of Torres et al., 2022 (Appendix A2.2). This filtering was only
used to identify the location of the submesoscale features and was not used when calculating the
Lagrangian trajectories.
Here, we define submesoscale features as locations with a Rossby number ≥ |0.75|. This
threshold was chosen based on a study by Molemaker et al., 2005 which demonstrated that unbalanced flow begins above Rossby numbers of 0.75 and subsequent work by McWilliams, 2016
and others showing that submesoscale features have Rossby numbers close to ±1. Three different
types of submesoscale features (cyclonic eddies, anticyclonic eddies, and fronts) were identified
based on the relative magnitude of Rossby number to normalized strain values (Torres et al., 2022)
(Figure 1c). The normalized local strain values were calculated as:
S =
r
du
dx −
dv
dy2
+
dv
dx +
du
dy2
f
(3.3)
Submesoscale fronts were defined as |ζ | ≥ 0.75 and S ≥ |ζ | (Torres et al., 2022). Submesoscale
eddies were defined as regions where |ζ | ≥ 0.75 but ¡S, where eddies with positive vorticity were
labeled cyclonic eddies and eddies with negative vorticity were labeled anticyclonic eddies. ζ
and S were calculated using u and v velocities from 10 m to avoid surface smoothing effects in
the vertical velocity. We quantified the fractional area occupied by submesoscale features at any
given timestep (spatial abundance of features) as the proportion of the regional surface area that
the three types of submesoscale features occupied. The number of submesoscale features at any
43
given timestep were also quantified using MATLAB’s bwcommcomp and regionprops functions
at each time step. The input for the functions was a binary array of the model domain where
ones corresponded to pixels where |ζ | ≥ 0.75 and zeros elsewhere. The bwcommcomp function
finds coherent features by identifying clusters of pixels with a value of one. Then, we used the
regionprops function to measure the area and feature pixel locations.
3.2.4 Spatial Analysis
To assess spatial patterns in Lagrangian particle dynamics, we broadly focus on four key regions:
North Nearshore, North Farshore, South Nearshore, South Farshore, with approximately equal
surface areas. We defined these regions based on distance from the coast, where the nearshore
regions extended 0-200 km and the farshore regions spanned 400-600km from the coast. The
North regions were between 38 - 41 and the southern regions spanned from 33.5 - 36.5 .
We quantified the encounter rate between the Lagrangian particles and submesoscale feature
within each broad region. A particle was assumed to be in a feature when the ζ and S calculated at
the latitude, longitude, and depth of the particle matched the criteria above (Figure 1d). We focus
our analyses on the mixed layer where the mixed layer was defined as the overall deepest point of
the mixing layer in a 24-hour period from the hourly output of the regional model. We calculated
the proportion of time spent in the mixed layer that a particle would encounter a submesoscale
feature as:
Particle-Feature Encounter Rate =
Σ time steps in a feature
Σ total time in the mixed layer
(3.4)
We tracked the total number of encounters between particles and submesoscale features, the
encounter frequency, and the type of submesoscale feature for each trajectory over the simulation
period. To further investigate the complex spatial dynamics governing particle-submesoscale feature interactions, we refined our analysis by subdividing the broader study regions into narrower
coastal proximity zones (0-25km, 25-50km, 50-100km, 100-200km, >200km). Subsequently, we
44
calculated the particle-feature encounter rate for each of these groups to examine potential variations in these interactions across different distances from the coast.
3.3 Results
To understand the physical environment experienced by particles in a submesoscale resolving
model and how this differs from the Eulerian fields, we advected over 2,000 Lagrangian particles in the CCS for almost 67 days (66.7) during the upwelling season. We assessed the vertical
and horizontal movements experienced by the particles. We compare the Lagrangian movements
and encounters with that of the overall model domain (Eulerian field). We present results from particles initialized at 100m, but findings were similar for particles released at 50m and 10m depths
(Appendix A2.1).
3.3.1 Eulerian Temperature and Mixing Variability
For the analyses presenting in this section, all Eulerian data was from the model output between
34.24 m and 39.725 m (mean of 37.0 m) which is approximately the mean depth of the mixed
layer experienced by the Lagrangian particles over the entire simulation time period . Overall, the
mean temperatures in the model were cool near the shore due to coastal upwelling and increased
offshore. Specifically, coastal temperatures in the north averaged 9.8 and in the south averaged
13.5 . Offshore temperatures increased to 12 in the northern part of the domain and 15 in the
south (Figure 2a).
Variability in model mixed layer temperature reflected the dynamic mesoscale flow of the region (Figure 2b), as expected. The temperature variance was organized into large-scale filaments
with low variability, and broad, windy swaths that outline eddy boundaries. Generally, the region
between 130 and 125 was the most variable region with up to ±- 2.4 (2σ) variability. Overall,
temperature variability in the Eulerian reference frame was an order of magnitude smaller than the
mean temperature values suggesting generally stable temperature conditions over the simulation
period.
45
Numerical models represent sub-grid scale advective processes as diffusive mixing (Kz). There
was a strong spatial gradient in mean diffusive mixing across the model domain (Figure 2c). Mixing was low near the coast, between 0.005 - 0.01 m2
/s, and increased offshore reaching a maximum of 0.035 m2
/s. The southern part of the model domain, near southern California, experienced
higher diffusive mixing, i.e., was more turbulent, than the northern region.
There was a distinct spatial pattern in the variability (2σ) of diffusive mixing rates. Locations
near the coast exhibited an order of magnitude less variability than offshore. Regions that were
marked with the highest mean diffusive mixing rates were the least variable and regions with
lower mean diffusive mixing rates were associated with the most variability in mixing (Figures 2d,
S3). This spatial pattern was especially pronounced near the coast where the standard deviation of
diffusive mixing was at least 7 times larger in magnitude than the mean diffusive mixing rate. Over
the whole model domain, the standard deviation was on average 2.8 ±2 (2σ) times higher than the
mean values.
3.3.2 Eulerian Submesoscale Feature Characteristics
The model represented three distinct categories of submesoscale features, fronts, cyclonic and
anticyclonic eddies. Across the model domain, we identified the occurrence of each feature type
and assessed the associated physical properties. To look at geographical differences, we divided
the study area into four quadrants: north/south and nearshore/farshore.
3.3.2.1 Submesoscale Feature Spatial Patterns
Overall, submesoscale features were rare covering on average 5% of the study area (Figure 3).
This is consistent with previous estimates of submesoscale prevalence (Liu et al., 2016; Ernst et
al., 2023). The nearshore regions had more features (combined 16% of regional areas) than the
farshore regions (combined 2.6% of regional area) – a 6.2 fold difference. The higher occurrence
of nearshore features is consistent with previous work which demonstrated that friction of flow
against the continental slope increases vorticity gradients, leading to more submesoscale feature
46
41
40
39
38
37
36
34
35
33
-132 -130 -128 -126 -124 -122 -120
0.035
0.03
0.025
0.02
0.015
0.005
0.01
0
41
40
39
38
37
36
34
35
33
-132 -130 -128 -126 -124 -122 -120
Kz [m2 ⁄ s] Kz mean/standard deviation
41
40
39
38
37
36
34
35
33
41
40
39
38
37
36
34
35
33
-132 -130 -128 -126 -124 -122 -120 -132 -130 -128 -126 -124 -122 -120
7
6
5
4
3
1
2
0
15
14
13
12
11
10
9
Temperature [oC] Temperature standard deviation [oC]
2.4
2.2
2.0
1.8
1.6
1.4
1.0
1.2
0.8
0.6
0.4
a b
c d
Figure 3.2: Model temperature and mixing conditions over the simulation period. a) Mean temperature at
39.72 5m depth for the period of March 27th -June 1st . b) standard deviation (2σ) at 39.725 m depth. c)
Mean diffusivity at 39.725 m for the whole simulation period. d) diffusivity standard deviation (2σ) divided
by the mean diffusivity.
47
generation (Lazar et al., 2018; Morvan et al., 2019; Dauhajre et al., 2017). We also observe an
increase in features in the southern region (combined 12% of total areas) compared to the northern
region (a combined 7% of total regional area) or 1.7 times more features in the south. The higher
occurrence of submesoscale features in the southern region is consistent with previous work which
showed that the southern portion of the CCS encompasses more eddy activity (Strub et al., 1991)
as a result of changes in wind stress curl (Dong et al., 2009; Winant et al., 1997) and topographical
changes (Dong et al., 2009; Caldeira et al., 2005; Dong et al., 2007) compared to the northern
region. This increased mesoscale eddy activity acts as a generation mechanism for submesoscale
features.
North
Onshore
North
Offshore
South
Onshore
South
Offshore
0
1000
2000
3000
4000
5000
6000
Total Feature Encounter by Lagrangian Particles
cyclonic
anticyclonic
front
0
0.02
0.04
0.06
0.08
0.1
0.12
Proportion of Region Area
cyclonic
anticyclonic
front
North
Onshore
North
Offshore
South
Onshore
South
Offshore
North
Onshore
North
Offshore
South
Onshore
South
Offshore
a b
Figure 3.3: Prevalence of submesoscale features. a) As a proportion of the surface area of each region b)
As a proportion of Lagrangian particle encounters per region.
Of the three types of submesoscale features, fronts were the most numerous in the CCS and
occurred almost twice as often (1.78 times) as cyclonic features and 2.23 times more than anticyclonic features. The prevalence of cyclonic vorticity and frontal features compared to anticyclonic
features was consistent with previous theoretical (Rudnick, 2001), observations (DiGiacomo et al.,
2001; Payandeh et al., 2023), and other modeling studies (Lin et al., 2020; Jing et al., 2020; Calil
et al., 2021). The similarities between submesoscale prevalence in the high-resolution model and
previous observational and modelling efforts provides confidence in subsequent analyses.
48
3.3.2.2 Submesoscale Feature Size
Different submesoscale features were associated with different spatial sizes in the model. Cyclonic
submesoscale eddies consistently had the largest spatial footprint, averaging 2.36 ±0.73 km2
(2σ),
followed by fronts (1.32 ±0.43 km2
), and anticyclonic eddies (1.15 ±0.42 km2
) (Figure A2.4).
The smaller size for anticyclonic eddies was consistent with regional observations (Payandeh et
al., 2023). Despite some overlap in size distributions, a two-sided t-test confirmed significant
differences in mean feature size among the different types (p<0.05). Additionally, for all feature
types, features in the southern regions were, on average, 29% larger (p<0.05) than their northern
region counter parts.
3.3.2.3 Submesoscale Velocities
The different submesoscale features exhibited distinct characteristic vertical velocities in both direction and magnitude (Figures 4 and S5). Anticyclonic features were associated with upwelling
due to the divergent nature of the features. Anticyclonic feature vertical velocity averaged 4.2
±10.5 m/day (2σ). Conversely, cyclonic features were associated with downwelling, consistent
with their convergent nature. Average cyclonic vertical velocities were -3.5 ±10.2 m/day (2σ),
where the negative value indicates downwelling. In general, fronts are characterized by weak
upwelling on the lighter side, and strong downwelling on the dense side (Capet et al., 2008b).
Because we were analyzing mean vertical velocities of the features, the result was net downward
vertical velocity for fronts. Submesoscale frontal features had the fastest average vertical velocities
-11.0 ±18.5 m/day (2σ). The background, (i.e. not within a submesoscale feature) average vertical
velocity was 0.1 ±9.9 m/day (2σ). The background vertical velocity was also significantly different from vertical velocities of the submesoscale features (two-sided ttest, p<0.05). The average
vertical velocities of the submesoscale features and the background did not vary significantly by
onshore/offshore or north/south regions.
49
-15 -10 -5 0 5 10
Eulerian Mean
-15
-10
-5
0
5
10
Lagrangian Mean
Average Vertical Velocities [m/day]
0.1 0.15 0.2 0.25 0.3
Eulerian Mean
0.1
0.15
0.2
0.25
0.3
Lagrangian Mean
Average Horizontal Velocities [m/s]
North Onshore
North Oshore
South Onshore
South Oshore Background
Cyclonic
Front
Anticyclonic
Figure 3.4: Mean vertical and horizontal velocities as calculated from the Eulerian reference frame compared the mean velocities encountered in the Lagrangian reference frame. The solid black line indicates the
Eulerian mean equal to the Lagrangian mean.
Horizontal velocities for all three feature types and the background varied less across the model
domain (Eulerian reference frame) than the vertical velocities, but all four values were still statistically different at the p<0.05 confidence level (two-sided ttest) (Figures 4 and S5). Anticyclonic
features were associated with the fastest horizontal velocities with an average of 0.17 m/s ±0.16
m/s (2σ) (Figure A2.6). Fronts were associated with the second fastest horizontal velocities with
an average of 0.14 m/s ±0.14 (2σ) and cyclonic features averaged 0.12 m/s ±0.12 m/s (2σ). Lastly,
the background magnitude of horizontal velocity was substantially less at 0.08 m/s ±0.09 m/s (2σ).
The background mean velocities were consistent with observational and model studies of the region (Matthews et al., 2009; Centurioni et al., 2008).
3.3.3 Lagrangian Trajectory Dynamics
The Lagrangian particles were advected in 3-dimensions throughout the model domain over the
66.7-day simulation. The most common advected path was in the southerly and offshore directions
50
Figure 3.5: Lagrangian particle trajectories over the entire simulation period. Trajectories are colored by
release latitude (see Figure 1 for release locations and final depths).
(Figure 5), as expected given the mean flow patterns in the region. On average, Lagrangian particles
traveled a total of 969 km ±289 km (2σ) over the 66.7 days. Overall, there was net upwelling of the
particles over the course of the simulation, also consistent with expectations based on large-scale
flow. Specifically, particles were released at 100 m depth and on average ended at 11 m ±22 m
(2σ).
The mixed layer shoaled over the simulation period transitioning from an average depth of 53.3
m to an average depth of 39.3 m, as expected for the springtime in this region (Figure 6). Here we
define the mixed layer as the deepest mixing point over a 24-hour period according to the K-Profile
Parameterization (Large et al., 1994). Within 10 days, 86% of the Lagrangian particles had reached
the mixed layer. Over the course of the simulation, the Lagrangian particles were within the mixed
layer 79% of the time. However, there were periods when the mixed layers in the region were
51
0 10 20 30 40 50 60 70
Simulation day
-200
-180
-160
-140
-120
-100
-80
-60
-40
-20
0
Mixed layer depth [m]
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Proportion of particles within the mixed layer
Mean Mixed Layer Depth
Mixed Layer Depth Range
Mean Particle Depth
Max Particle Depth
Prop. Particles in
mixed layer
Figure 3.6: Mean mixed layer depth across all Lagrangian particles over time (blue). Blue shading represents the range of the mixed layer depth across all particles at each timestep. Also shown is the mean depth
of the Lagrangian particles over time (black) and the maximum depth across all particles at each time step
(dotted black line). The green line represents the proportion of particles that are within the mixed layer.
Changes in the proportion of particles in the mixed layer results from mixed layer dynamics opposed to
vertical advection of the particles.
deep and over 90% of the Lagrangian particles are within the mixed layer (e.g. simulation days
17 and 40) and other times when the mixed layers were shallow and very few particles were found
within the mixed layer (e.g., day 27 when only 14% were in the mixed layer). The variability in
the proportion of particles in the mixed layer was primarily driven by variability in the mixed layer
depth rather than vertical advection of the particles.
52
3.3.4 Reference Frame Impact on Temperatures Experienced
The mean environmental conditions experienced by the Lagrangian particles was significantly
different than the variability quantified in the Eulerian reference frame – as might be expected
(two-sided ttest, p<0.05) (Figure 7a). On average, Lagrangian particles experienced lower mean
temperatures over the 66 days than the corresponding Eulerian values at the final location of the
particle (range -2.7 - 1.5 ). This is consistent with the general offshore movement of the trajectories and the cooler temperatures nearshore due to upwelling (Figure 2). Similarly, particles that
traveled more than 400 km offshore all experienced lower average temperatures than the corresponding Eulerian mean and particles that traveled more than 200 km in the onshore direction had
higher Lagrangian means than the Eulerian fixed point (Figure 7a).
The Lagrangian particles also experienced different variability (2σ) in temperature relative
to the Eulerian reference frame (two-sided ttest, p<0.05) (Figure 7a). On average, Lagrangian
particles experienced 2.4 times (1 ) more variability than the corresponding Eulerian values (range
-0.5 - 1.9 ). Similar to the differences in mean temperatures experienced, variability in temperature
varied with distance traveled relative to the coast. Generally, as particles traveled offshore, the
Lagrangian temperature variability increased relative to the Eulerian variability.
3.3.5 Reference Frame Impact on Mixing Rates Encountered
Previous work has suggested that mixing rates in the Eulerian reference frame are larger than those
encountered in the Lagrangian reference frame (Baudry et al., 2018; Paparella et al., 2020; Hellweger et al., 2007). However, for our model domain, there was not a significant difference between
the mixing experienced in the two reference frames with 53% of the Lagrangian particles experienced higher mean mixing rates than the Eulerian reference frame and 47% experienced lower
mean mixing rates (Figure 7b, -0.02 - 0.01 m2
/s). However, a spatial relationship was observed.
Particles that traveled greater distances offshore tended to have lower and less variable mean mixing rates than the Eulerian reference frame while particles advected onshore had higher and more
variable mixing rates.
53
-1000 -800 -600 -400 -200 0 200 400
-0.025
-0.02
-0.015
-0.01
-0.005
0
0.005
0.01
0.015
-0.5
0
0.5
1
Kz [m2⁄s]
Lagrangian mean - Eulerian mean
(Lagrangian
σ - Eulerian
σ)/Eulerian
σ (Lagrangian
σ - Eulerian
σ)/Eulerian
σ
-1000 -800 -600 -400 -200 0 200 400
ΔDistance from coast over simulation [km]
-3
-2.5
-2
-1.5
-1
-0.5
0
0.5
1
1.5
Lagrangian mean - Eulerian mean
-0.5
0
0.5
1
1.5
2
2.5
3
oshore advection onshore advection
Temperature [oC]
ΔDistance from coast over simulation [km]
oshore advection onshore advection
Figure 3.7: Difference between Lagrangian and Eulerian reference frames as a function of particle advection. (a) As Lagrangian particles are advected offshore, the mean temperature along the trajectory decreases
relative to the Eulerian mean. However, the temperature variability experienced by the Lagrangian particles
advected offshore is higher than the corresponding Eulerian variability. (b) Same as panel (a) but showing
differences in diffusive mixing values between the two reference frames. Although the difference in means
between the two reference frames follows the same spatial pattern as temperature, the Lagrangian variability
decreases relative to the Eulerian reference frame with offshore advection.
54
3.3.6 Reference Frame Impact on Velocity
Our findings indicate that submesoscale features, despite their intense vertical velocities, had a
limited impact on the vertical advection of Lagrangian particles in the CCS. This was primarily
due to the brief residence times of particles within these features.
To understand the factors influencing these short residence times, we compared Eulerian and
Lagrangian vertical velocities. Specifically, we tested if the vertical velocities encountered by the
Lagrangian particles were larger in magnitude due a sampling bias in the Lagrangian reference
frame which would result in even shorter residence times than expected based on Eulerian means.
In the Eulerian reference frame, we observed a clear pattern in the vertical velocities associated
with different feature types: background velocities were slightly upwelling, anticyclonic features
were predominantly upwelling, and cyclonic features and fronts were both downwelling, with
fronts exhibiting stronger downwelling (Figure 4).
However, in the Lagrangian reference frame, this pattern was less distinct, with more noise and
overlap in the vertical velocities associated with each feature type. While Lagrangian particle advection showed a wider range of extreme velocities, there was no systematic bias towards stronger
vertical motion within any feature type relative to the Eulerian means (Figures 4 and S5).
These results suggest that the short residence times of particles within submesoscale features
are primarily driven by the strong vertical velocities associated with these features as also seen
in the Eulerian reference frame, rather than any systematic sampling bias in the Lagrangian reference frame. Furthermore, the relationship between feature type and vertical advection direction
observed in the Eulerian reference frame is not as clear in the Lagrangian reference frame.
While most studies on submesoscale advection focus solely on the vertical component, the
strong surface divergence often associated with these features (McWilliams, 2016) suggests a potential role in horizontal advection as well. Our analysis confirms this, revealing significant differences in horizontal velocity magnitudes between different feature types and the background flow.
In the Eulerian reference frame, anticyclonic features exhibited the fastest horizontal velocities.
However, the Lagrangian reference frame showed that particles experienced the highest velocities
55
2 4 6 8 10 12
Mean Time in Feature [hr]
0
0.2
0.4
0.6
0.8
1
1.2
Kernel Density Estimate
Cyclonic Features
North Nearshore
North Farshore
South Nearshore
South Farshore
2 4 6 8 10 12 0
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
0.18 Anticyclonic Features
2 4 6 8 10 12 0
0.2
0.4
0.6
0.8
1.0
1.2
1.4
1.6 Fronts
North Nearshore
North Farshore
South Nearshore
South Farshore
Mean Time in Feature [hr]
Kernel Density Estimate
Mean Time in Feature [hr]
North Nearshore
North Farshore
South Nearshore
South Farshore
Kernel Density Estimate
Figure 3.8: Kernel density estimates (analogous to probability distributions) of the time spent in a submesoscale feature by the Lagrangian particles. The data are separated by region and data type. Generally,
particles spent the longest time in cyclonic features. However, the mean time spent in a feature is 1.6 hours,
much shorter than the lifetime of the features.
when interacting with submesoscale fronts, averaging 0.24 m/s ±0.31 m/s (2σ) - over 150% faster
than the Eulerian average for fronts (Figure A2.6). In fact, all Lagrangian horizontal velocities
exceeded their Eulerian counterparts (Figure 4).
Notably, particles within cyclonic features also experienced significantly faster horizontal velocities (0.22 m/s ±0.30 m/s) than the Eulerian average. While anticyclonic and cyclonic features
showed similar mean velocities in the Lagrangian frame, they remained statistically distinct at the
95% confidence interval. As seen with the vertical velocities, there was a much clearer separation of horizontal velocities by feature type in the Eulerian reference frame than the Lagrangian
reference frame (Figure A2.6).
Consistent with the Eulerian data, background flow velocities in the Lagrangian frame were
slower than within features, yet still notably faster than their Eulerian counterparts (0.17 m/s ±0.21
m/s), highlighting the dynamic nature of particle advection in this region.
The short residence times of Lagrangian particles within submesoscale features were primarily
a consequence of the intense vertical and horizontal velocities associated with these features. This
dynamic advection, while limiting the direct impact of submesoscale features on vertical transport,
highlights their significant role in shaping the overall Lagrangian flow field in the CCS.
56
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09
Occurrence Rate by Total Area
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
Encounter Rate by drifters in mixed layer
0
50
100
150
200
250
Distance from coast [km]
cyc-north
anticyc-north
front-north
cyc-south
anticyc-south
front-south
Figure 3.9: Comparison of Eulerian occurrence rate by area to the Lagrangian encounter rate by time spent
in the mixed layer.
57
3.3.7 Lagrangian Encounters by Submesoscale Feature Type
We quantified the total number of submesoscale features encounters by the Lagrangian particles,
the number of timesteps each particle spent in each feature, and the frequency with which the
Lagrangian particles encountered submesoscale features. An encounter was defined as points along
the trajectory when the absolute value of the Rossby number at the particle’s location was ≥ 0.75.
We also classified each encounter based on the type of submesoscale feature (i.e., cyclonic eddy,
anticyclonic eddy, or front) based on the normalized strain rate at the particle’s location. In general,
Lagrangian trajectories encountered more submesoscale features when they were closer to shore
and fewer features when further from shore, consistent with the Eulerian reference frame (Figure
3). Similarly, cyclonic and frontal features occupied the greatest fractional spatial area in the
Eulerian analysis and were the most encountered type of features by the Lagrangian particles.
Although submesoscale features can last on the order of days, Lagrangian particles were directly in a submesoscale feature for on average only 1.6 hrs ±2.2 hrs (2σ) at a time with a range
of 1-15 hours (Figure 8). Overall, we observed that Lagrangian particles were advected out of features on timescales shorter than typical lifespans of submesoscale features. For the majority of the
model domain (all regions except the north farshore), Lagrangian particles stayed within cyclonic
features for the longest period of time of 1.7 hrs ±2.5 hrs (2σ), which is consistent with cyclonic
features having the slowest mean vertical velocity. In the north farshore region, Lagrangian particles spent the most time in anticyclonic features with an average of 2.0 ±4.4 hrs (2σ) in each
feature. For all regions, trajectories were within fronts for the shortest amount of time averaging
just 1.4 hrs ±1.5 hrs (2σ). Fronts were associated with the fastest downward velocities (Figure
A2.5) suggesting that these short residence times within fronts were due to vertical advection out
of the feature.
If encounter and residence time within submesoscale features were stochastic, we would assume that the time spent by Lagrangian trajectories in submesoscale features would match the
occurrence rate (fractional spatial area occupied) of the features. However, if trajectories became
trapped within features, we would expect an apparent over sampling of the features. Conversely,
58
if vertical or horizontal advection within the features pushed the Lagrangian particles out of the
features (e.g. subduction below the front) then we would expect the Lagrangian trajectories would
under sample features relative to their occurrence. Given the above analysis of particle residence
time within features, we expected to find that the trajectories under sampled submesoscale features.
We observed significant under sampling, particularly for particles ¿150km from the coast. A statistically significant, linear relationship between Lagrangian encounter rates and feature occurrence
rates (slope = 0.1889, R2 = 0.83, p <0.05) indicated that offshore particles encountered features
80% less frequently than their occurrence. This under-sampling is consistent with stronger vertical
velocities for offshore submesoscale features (Figures S7 and S8) which act to decrease the residence time of trajectories within features relative to the nearshore (Figure 8). In contrast, nearshore
particles, exhibited encounter rates commensurate with feature occurrence, likely due to extended
residence times within these features relative to their offshore counterparts (Figure 8). This spatial variability in encounter rates highlights the complex interplay between feature distribution and
particle residence time in different oceanic regimes.
Our analysis indicated that submesoscale features only occupy a small fraction of the spatial area and are encountered by Lagrangian trajectories, on average, less than 5% of the time.
However, previous studies have suggested that the rapid vertical velocities associated with submesoscale features might play an important role in biogeochemical cycles even if the prevalence of
these features is low (Liu et al., 2016). We assessed this by quantifying the relative contribution
of submesoscale vertical advection to the total vertical advection experienced by the particle. The
contribution of submesoscale vertical advection to the total vertical advection was less than 5%,
consistent with the Lagrangian encounter rates and time spent within the features. This indicates
that the influence of submesoscale dynamics on the vertical movements of our trajectories, while
present, was not disproportionate to their spatial and temporal prevalence.
59
0 100 200 300 400 500 600 700
Drifter Number
0
0.01
0.02
0.03
0.04
0.05
Combined Vertical Advection Proportion
cyclonic
front
anticyclonic
Figure 3.10: Proportion of total vertical advection resulting solely from advection due to submesoscale
features. The analysis was performed across all Lagrangian particles. Lagrangian trajectories with submesoscale proportions of less than 1.9*10−6
are shown.
3.4 Discussion
We assessed differences in environmental variability in the Lagrangian reference frame compared
to the Eulerian reference frame using 3D Lagrangian trajectories within a high-resolution (500 m)
physical ocean model the of the CCS region. We found that not only are the mean states different
between the two reference frames, but also the variability is different. Specifically, for temperature and diffusive mixing rates, there is a strong spatial component related to the onshore/offshore
advection of the particle that controls the difference between the reference frames. Our results
challenge the view that the Eulerian reference frame consistently overestimates Lagrangian mixing rates. Rather, our results suggest that the two reference frames are significantly different but
that approximately half the time the Eulerian reference frame has higher mixing rates than the Lagrangian and vice versa. Whether the Eulerian reference frame over- or underestimates Lagrangian
mixing was a function of horizontal advection. When Lagrangian particles moved from a region
of low mixing rates (e.g. near the coast) to regions of higher mixing rates (e.g. offshore), the
Lagrangian mean mixing rates were lower than the corresponding Eulerian mean and vice versa.
The environmental differences between reference frames have biogeochemical implications.
Typically, biogeochemical models are run in the Eulerian reference frame whereas phytoplankton
60
experience the environment in the Lagrangian reference frame. Because the mean temperature
and mixing rates between the two reference frames differ, under constant, mean conditions alone,
phytoplankton growth dynamics are likely to be different between Eulerian and Lagrangian based
biogeochemical models (e.g. Zaiss et al., 2021). Taking environmental variability into account will
result in overall less growth for both reference frames due to non-linear growth dynamics (Jensen,
1906). The differences in phytoplankton growth dynamics between Eulerian and Lagrangian biogeochemical models are going to be even larger when the difference in reference frame variability
is considered. However, Lagrangian particle simulations are time and computer power intensive
so biogeochemical models, especially large-scale models, operate in the Eulerian reference frame.
Our results highlight the need for future work to relate Eulerian and Lagrangian environmental
conditions to one another so that Eulerian biogeochemical models can better capture true phytoplankton dynamics.
One way in which establishing Eulerian-Lagrangian environmental relationship may be useful
is when trying to incorporate the impact of fine-scale physical dynamics of the CCS into coarser
grid, large-scale models that do not resolve submesoscale features. Eulerian estimations suggest
that submesoscale features cover 5% or less of the coarse-grid cell. When accounting for the
changes in vertical and horizontal velocity within that grid cell, a 5% adjustment seems logical.
However, our results indicate that Lagrangian particle interactions with submesoscale features differs from what would be predicted in the Eulerian reference frame. An Eulerian based estimate
would underestimate Lagrangian particle-feature encounter rates nearshore and over-estimate the
encounters offshore. One reason why this matters is because observation-based work has found
that in the CCS, it takes 3-7 days of surface residence time to start accumulating biomass (Dugdale
et al., 2006; Dugdale et al., 1989; Wilkerson et al., 2006). Here, we show that Lagrangian particles
are more likely to encounter downwelling features (fronts and cyclonic features) than upwelling
features (anticyclonic features). From submesoscale vertical advection alone, Lagrangian particles are more likely to be subducted out of the surface waters preventing the particle reaching the
required residence time for biomass accumulation to begin. The frequency in which Lagrangian
61
particles encounter submesoscale features will influence whether biomass can accumulate along
the trajectory. Specifically, encounters with downwelling features more frequently than once every
3-7 days will prevent biomass from accumulating along that trajectory. In the nearshore region, the
impact of submesoscale subduction would be underestimated which would permit more biomass
to accumulate in the Eulerian reference frame than would likely occur in the Lagrangian reference
frame and vice versa in the offshore regions.
3.5 Conclusions
Our study highlights the impact of reference frame on environmental mean conditions and variability experienced. Additionally, we show that trying to account for Lagrangian conditions in
the Eulerian reference frame is not a one-to-one approximation and varies spatially, specifically
with distance from the coast. As such, incorporating Lagrangian environmental dynamics within
Eulerian-based biogeochemical models is crucial and complex. The results presented here will
hopefully act as a guide for future research to consider the Lagrangian experience in Eulerian
models. Establishing a relationship between Lagrangian and Eulerian environmental conditions in
other seasons is crucial for developing a comprehensive framework for combining the two reference frames.
In addition to differences in simulated reference frames, we need more observational data about
the environmental differences between Eulerian and Lagrangian perspectives. Specifically, there
needs to be more truly Lagrangian observations to accompany Eulerian observations. Autonomous
vehicles and gliders offer some insight to Lagrangian environment variability but are not completely passive observational techniques. Establishing a relationship between what is measured
in an Eulerian reference frame and what is encountered along Lagrangian pathways is critical to
include the impact of fine-scale features in biogeochemical models, particularly those used to investigate changes in the global carbon cycle under warming conditions.
62
Chapter 4
A Lagrangian Lens on Biophysical Interactions in the
California Current System
Abstract
This study presents a novel approach for investigating the biophysical controls on phytoplankton
biomass accumulation. We developed a one-dimensional, Lagrangian, water column model and
coupled it with a high-resolution (500m) submesoscale-resolving physical ocean model. With this
model framework, we quantified the balance between physical driven biomass losses and gains,
and biomass changes that were biologically driven through growth and mortality. This study reveals the complex interplay between mixed layer depth, photic zone dynamics, and particle depth
as they impact biomass accumulation. We demonstrate that the highest phytoplankton growth
rates do not always translate to the most biomass accumulation. Rather, diffusive mixing plays a
key role by decreasing surface biomass concentration through dilution which exports biomass out
of the surface. Our work highlights the importance of studying phytoplankton dynamics in the
Lagrangian reference frame.
4.1 Introduction
Phytoplankton are crucial to the global carbon cycle. Half of all global photosynthesis is performed
by phytoplankton which acts to draw down atmospheric carbon dioxide. The rate of photosynthesis and resulting accumulation of biomass is a function of the environmental history experienced
by the phytoplankton (Hellweger et al., 2016; Baudry et al., 2018; Boyd et al., 2016; Doblin et al.,
63
2016; Zaiss et al., 2021). As phytoplankton are largely passive drifters, the environmental variability they experience occurs in the Lagrangian reference frame, along the water mass trajectory.
However, typically biogeochemical models simulate phytoplankton dynamics in the fixed point,
Eulerian reference frame. Previous work has shown that the history of variability encountered
by a Lagrangian path differs from the Eulerian variability experienced at the same location (e.g.
Doblin et al., 2016). For instance, tracking temperature along Lagrangian trajectories results in
average temperatures for the moving parcel that can be quite substantially different from the average temperature at the final location. In addition, there are also differences in the environmental
variability experienced along a Lagrangian trajectory compared to the end location. Differences
in environmental variability (the duration and frequency of changes in environmental conditions)
have been shown to impact phytoplankton growth and biomass accumulation rates (Zaiss et al.,
2021; Baudry et al., 2018; Bernhardt et al., 2018; Schaum et al., 2018). To fully understand the
physical and biological mechanisms that impact phytoplankton growth, it is imperative to capture
variability in the Lagrangian reference frame.
There has recently been an increase in studies focused on phytoplankton growth dynamics from
the Lagrangian perspective using satellite and float data, as well as model simulations (Jonsson et ¨
al., 2016; Lehahn et al., 2018; McKee et al., 2023; Frischknecht et al., 2018; Zaiss et al., 2021;
Baudry et al., 2018). These methods have made great strides in understanding the mesoscale
(>100km, months) processes that impact phytoplankton accumulation rates. Specifically, it has
been shown that mesoscale strain can relieve light limitation through mixed layer shoaling and
promote phytoplankton growth (McKee et al., 2023). Additionally, the rate and magnitude of
temperature variability can play a significant role in phytoplankton community growth rates (Zaiss
et al., 2021).
However, fine-scale physics, generally defined here as smaller than mesoscale (submesoscale,
O1-10km, days), are increasingly recognized for their role in nutrient transport and carbon export
(Freilich et al., 2019; Brannigan, 2016; Cao et al., 2024; Chrysagi et al., 2021; Liu et al., 2018;
Taylor et al., 2020). Submesoscale eddies and fronts may be particularly relevant for phytoplankton
64
growth dynamics due to their rapid vertical velocities, up to 100m/day (Mahadevan et al., 2006).
These rapid vertical velocities can act to inject nutrients into the euphotic zone, upwell light limited
phytoplankton, or subduct biomass out of the surface waters which also acts to increase light
limitation for the subducted phytoplankton (Levy et al., 2012b; Mahadevan, 2016). ´
To understand the impact of fine-scale physics on phytoplankton growth and biomass accumulation, modeling studies have conventionally relied on comparing coarse, non-submesoscale
permitting simulations to output from high-resolution, submesoscale permitting or resolving models. However, increasing model resolution changes flow patterns and better resolves larger-scale
features which complicates direct comparison between models (Balwada et al., 2021; Busecke et
al., 2019; Levy et al., 2010). This method of comparison also neglects the potential differences ´
in phytoplankton growth rates caused by differences in Eulerian versus Lagrangian environmental
variability. Thus, we still lack fundamental understanding of how fine-scale physical dynamics
impact phytoplankton growth and biomass accumulation in large part because we have not had the
proper tools available to study bio-physical interactions at the correct temporal and spatial scales.
In addition to submesoscale features, which have been the focus of most of the primary research, there are other fine-scale physical dynamics that are resolved in high-resolution models
(Trotta et al., 2017). One such dynamic is the significant (up to an order of magnitude) increase in
vertical diffusivity, or mixing rate, that these models capture compared to coarser resolution models (Liu et al., 2021). Diffusive mixing plays a key role in shaping the environment experienced
by phytoplankton by altering key factors such as nutrient availability and light exposure. For example, (Baudry et al., 2018) demonstrated that stratified or intermediate levels of mixing leads to
Eulerian-based phytoplankton models overestimating phytoplankton biomass compared to the Lagrangian reference frame. Understanding the coupled influence of vertical advection and mixing
rates in shaping phytoplankton environmental history and, consequently, biomass accumulation,
is key to fully grasping the impact of fine-scale bio-physical interactions. Additionally, internal
gravity waves, which are resolved better in high-resolution models as well, have similar characteristic vertical velocities and spatiotemporal scales as submesoscale features. The vertical advection
65
associated with internal gravity waves can contribute to changes in light limitation and impact phytoplankton growth (Evans et al., 2008; Garwood et al., 2020). Further, breaking internal gravity
waves can result in upward mixing of nutrients to the sunlit surface waters (Sandstrom et al., 1984;
Lucas et al., 2011a). Therefore, distinguishing the relative contributions of all fine-scale physical
processes is imperative for a full understanding of biophysical controls on phytoplankton biomass.
The California Current System (CCS) is a great system for studying bio-physical interactions.
The CCS is characterized by wind-driven upwelling in the spring which creates strong spatial gradients in temperature, salinity, nutrients (Abbott et al., 1990; Castelao et al., 2006; Chavez et al.,
1991; Checkley et al., 2009; Deutsch et al., 2021; Lynn et al., 1987) and ultimately phytoplankton
growth (Taylor et al., 2012; Taylor et al., 2018; Lynn et al., 2003; Hood et al., 1990; Dugdale et al.,
1989). In addition, mesoscale eddies formed by a north-south baroclinic jet results in the formation
of submesoscale features (Capet et al., 2008a; Payandeh et al., 2023). In addition to rapid vertical velocities due to submesoscale features, internal gravity waves contribute significantly to the
physical dynamics in the CCS (Torres et al., 2018; Torres et al., 2022). Therefore, characterizing
the bio-physical interactions with submesoscale features in the CCS is imperative to furthering our
understanding of the roles that fine-scale physics plays in phytoplankton biology.
4.2 Methods
To investigate the impacts of submesoscale dynamics on phytoplankton growth and biomass accumulation, we developed a novel modeling framework which couples Lagrangian particle trajectories derived from a high-resolution, submesoscale-permitting ocean model (Torres et al., 2022)
with a phytoplankton biology model.
4.2.1 Physical Model
Here, we use a high resolution, physical oceanographic, regional model simulation for the CCS.
See Chapter 2 Methods for model details.
66
4.2.2 Lagrangian Trajectories
We released at total of 2,065 Lagrangian particles (Figure 1A) to advect over the entire simulation
period (Figure 1B). See Chapter 2 Methods for details on particle initialization and advection.
After filtering trajectories that went outside the model domain, there were 1,463 67-day particle
trajectories used in the final analyses.
C D
A B
15
50
100
150
200
Figure 4.1: : Examples of Lagrangian drifter trajectories. A) Starting locations of all the Lagrangian
drifters at the first time step. drifters were initialized at a depth of 100m and are colored according to
starting latitude. The background is colored by sea surface temperature [] at the first time step. B) Final
location of the drifters on day 66. The color represents the latitude that the drifters started at in (A). The
size of the drifter represents the depth of the drifter. Larger point size mean deeper drifters as specified in
the legend (meters). C) Trajectory of a Lagrangian drifter through time. The background is colored by the
relative vorticity at day 54.5 and the trajectory is colored by simulation day. At simulation day 54.5, the
Lagrangian drifter encounters a cyclonic frontal feature, marked by a star. (D) The depth through time of
the trajectory in (C). The line is colored by the relative vorticity that the Lagrangian drifter encounters. The
initial encounter with the frontal feature is marked by a start at simulation day 54.5. The drifter is rapidly
subducted by the frontal flow.
67
4.2.3 Phytoplankton Model
To simulate phytoplankton growth and response to fine-scale environmental variability, we implemented a modified cellular quota model. This model, based on Ward et al., 2017 with some
modifications, simulates phytoplankton chlorophyll (Chl) concentration and biomass (C). Phytoplankton growth is limited by light, nutrients, and temperature. Chlorophyll content varies with
light levels and is tied to nitrate uptake. The flexible Chl to C ratios of the quota model enables us
to distinguish whether phytoplankton growth is primarily limited by light due to vertical advection
or nutrient availability. A table of all variables and their units can be found in Tables A3.1-A3.3 in
the appendix.
4.2.3.1 Biomass and Nitrogen Equations
The Lagrangian water-column model solved for water column nitrate and ammonia concentrations
and photosynthetically active radiation (PAR) with depth. Phytoplankton growth and mortality
was calculated based on these external variables and the change in three separate internal biomass
pools over times were tracked: phytoplankton carbon content, phytoplankton Chl content, and
phytoplankton nitrogen content. The model also tracked the total number of cells (X).
The phytoplankton carbon biomass (C) in each model grid box was calculated as
dC
dt
+w
dC
dt
+Kz
d
2C
dz2
= µ ∗C − MortC (4.1)
The second two terms on the left represent the impact of physics on phytoplankton biomass
where w
dC
dt is the advective loss of biomass and Kz
d
2C
dz2
is the impact of diffusive mixing on biomass.
The right-hand side terms represent growth and mortality dynamics where µ is the daily growth
rate and MortC is mortality due to grazing and cell death. The number of phytoplankton cells (X)
was then calculated based on the amount of phytoplankton carbon in each grid box using a fixed
relationship for carbon content per cell ratio (6.04e-5 mmol C/cell).
dChl
dt
+w
dChl
dt
+Kz
d
2Chl
dz2
= Chlsyn ∗X − MortC ∗ChlC (4.2)
68
38
37.5
36
36.5
35.5
35
34.5
34
-128.5 -128 -127.5 -127 -126.5 -126 -125.5
0
-40
-80
-120
-160
-200
0 10 20 30 40 50 60
0
4
8
12
16
20
Carbon [mmol/m3]
Depth [m]
Latitude
Longitude
Simulation Day
Day 0
Day 67
Day 0
Day 67
0
-40
-80
-120
-160
-200
Depth [m]
Figure 4.2: : Example of Lagrangian phytoplankton model. (Top) Idealized latitude and longitude of
a Lagrangian particle over the whole simulation. (Middle) Depth profile through time of the idealized
trajectory in the top panel. (Bottom) 1-D water column model with phytoplankton biology advected along
the Lagrangian particle trajectory over the entire simulation.
69
where Chlsyn is the Chl synthesis rate per cell. . The Chl loss term is the carbon loss rate
multiplied by the Chl:C ratio. The change in N biomass over time was calculated as:
dN
dt
+w
dN
dt
+Kz
d
2Nl
dz2
= (VNO3 +VNH4) ∗X − MortC ∗NC (4.3)
where VNO3 and VNH4 are the cellular uptake rates of nitrate and ammonia, respectively. The
nitrogen biomass loss term is the cellular N:C multiplied by the carbon mortality. Nitrate concentrations in the biology model only changed due to cellular uptake and physics and were calculated
as:
dNO3
dt
+w
dNO3
dt
+Kz
d
2NO3
dz2
= −VNO3 ∗X (4.4)
Ammonia concentrations changed as a function of physics, uptake, and remineralization of
nitrogen biomass mortality calculated as:
dNH4
dt
+w
dNH4
dt
+Kz
d
2NH4
dz2
= (0.1 ∗ MortC ∗NC −VNH4) ∗X) (4.5)
4.2.3.2 Cellular C Quota
The cellular carbon quota (Qc [mmol C/cell]) was calculated as a function of cell volume (Vcell
[µm
3
]) as follows:
Qc = 3.6 ∗ 10−11 ∗V
0.94
cell (4.6)
4.2.3.3 Nitrogen Uptake Rate
Nitrogen is considered the primary limiting nutrient in the model. This is consistent with other
high-resolution phytoplankton models of the CCS (Kessouri et al., 2020; Damien et al., 2023). Our
model simulates both new primary productivity (from nitrate) and regenerated (from ammonia).
Nitrate nitrate uptake rates (mmol N/cell/day) were calculated as:
70
VNO3 = Vmax,NO3 ∗ γ
T
∗ kinhib (4.7)
where kinhib is the ammonia inhibition, = Vmax,NO3 is the maximum rate of nitrate uptake per
cell, and γ
T
is the temperature limitation factor. Following (Ward et al., 2016) the maximum uptake
rate scales with size as:
Vmax,NO3 = Qc ∗ 0.51 ∗V
−0.27
cell (4.8)
If ammonia was available, it was the preferred nutrient before nitrate. Therefore, kinhib is
calculated as:
kinhib = exp(−ψ ∗NH4) (4.9)
where ψ is an ammonia inhibition constant and NH4 is the concentration of ammonia in
mmol/m3
.
The ammonia uptake in mmol NH4 per cell per day is a function of temperature limitation and
ammonia availability:
VNH4 = Vmax,NH4 ∗ γ
T
∗
NH4
(NH4 +kNH4
(4.10)
where NH4 is the concentration of ammonia in mmol/m3
, Vmax,NH4 is the maximum uptake rate
of ammonia in mmol N/cell/day and scales with cell size as:
Vmax,NH4 = Qc ∗ 0.26 ∗V
−0.27
cell (4.11)
γ
T
is the temperature limitation calculated as:
γ
T = τ ∗ exp
A∗
1
T
−
1
To
(4.12)
71
where τ is a unitless constant, A is a constant with units of Kelvin, To is the maximum growth
temperature in Kelvin, and T is the temperature at the current time step.
4.2.3.4 Photosynthesis and Respiration
The photosynthesis rate is modified from a maximum photosynthetic rate altered by light limitation. Photosynthesis rate was calculated as
P = Pmax ∗ γ
L
(4.13)
where the maximum photosynthesis rate, Pmax, is calculated as
Pmax = µmax ∗ γ
T
∗ γ
N
(4.14)
where µmax is the phytoplankton maximum growth rate γ
N is the nutrient limitation that follows
Michalis-Menton kinetics:
γ
N =
NO3
(NO3 +kNO3
(4.15)
γ
L
is the limitation due to light. The light limitation is a function of cellular chlorophyll calculated following a modified version of Li et al., 2010 as
γ
L =
1−exp
−α ∗ChlC∗ I
Pmax
∗ exp
−β ∗ I
Pmax
(4.16)
where α is the initial photosynthesis-irradiance slope, β is photoinhibition parameter, Chl : C
is the cellular chlorophyll to carbon ration, and I is the irradiance. γL was then used to calculate
photosynthesis rate (P, eq 13). Additionally, at each time step for each trajectory, light limitation
values (γL) were normalized so that maximum value (least limited) within the water column was
1. This value was typically located at the surface but due to photoinhibition in the model, could
sometimes be slightly below the surface level.
Net phytoplankton growth was then calculated as:
72
µ = P−Rb (4.17)
where Rb is the basal respiration cost which is a constant.
For this model we used a linear, temperature dependent mortality rate (MortC) which was
calculated as:
MortC = 0.001 ∗ γ
T
∗C (4.18)
Other forms of mortality including quadratic were tested in model development. Linear mortality was chosen because it provided the best fit to the observed distributions of carbon, chlorophyll,
and nitrate.
4.2.3.5 Chlorophyll Synthesis and Degradation
Chlorophyll was synthesized in the model as a function of nitrate uptake following Geider et al.,
1996 as:
Chlsyn = ChlNmax ∗
P
α ∗ChlC∗ I +β ∗ I
∗ (VNO3 +VNH4) (4.19)
where ChlNmax is the maximum cellular chlorophyll to nitrogen ratio.
4.2.4 1-D water column model
Along each of the Lagrangian particle trajectories, we simulated a one-dimensional water column
from the surface down to 200m using 51 evenly spaced grid boxes (dz=4m) (Figure A3.1). At each
time step, the advection, diffusion, and temperature from the surface to 200m were extracted from
the physical model and linearly interpolated to align with the vertical grid spacing of our water
column model. These values were then used to force the biogeochemical model and simulate
substrate changes due to advection and diffusion, represented in the left-hand side of equations
1-5.
73
The biogeochemical model was integrated forward in time, driven by physical variables derived
from the CCS regional model. To ensure numerical stability, grid cells with diffusivities exceeding
0.026 m2
/s were instantaneously mixed, resulting in homogeneous model state variables (concentrations of nutrients, Chl, and biomass) over the region with elevated rates of diffusive mixing.
The boundary conditions conditions for the model was set as follows. For nitrate, the concentration for bottom grid box was set using the established linear relationship between temperature
and nitrate for the region(Lucas et al., 2011b). The bottom boundary concentrations for phytoplankton carbon, nitrogen and chlorophyll were set to minimum values of 0.1 mmol C/m3
, 0.02
mmol N/m3
, and 0.02 mg Chl/mmol C (at the minimum Chl:C).
Light levels at the surface (Io) were calculated based on time of day and latitude following
(Brock, 1981). First the declination was calculated from the number day N of the year as
D = 23.45 ∗ sin360 ∗ (284+N)./365. (4.20)
Then the hour-angle, in degrees, was calculated based on the hours since midnight (T) as
W = (T −12) ∗ 15 (4.21)
Now using the declination, the drifter latitude L in degrees, and the hour-angle, the cosine of
the Zenith angle, in degrees, was calculated as:
cosZ = sinD∗ sinL+cosD∗ cosL ∗ cosW. (4.22)
Finally, the clear sky radiation, PARo (W/m
2
), was computed as
PARo = (1121 ∗ (cosZ −0.08251)). (4.23)
Light attenuation with depth was derived as
74
I = Io ∗ exp(−kd ∗ dz) (4.24)
where Io is the surface light, dz is the depth of the grid boxes, and kd is the attenuation due to
seawater and chlorophyll. The attenuation was calculated as
kd = ksw+kchl ∗Chl (4.25)
where ksw and Kchl are constants and Chl is the chlorophyll concentration. The kd was calculated for each depth as a function of the light and Chl in the layer directly above.
Initial conditions for the water column model were set following a two-step, spin up process.
The first step initialized the water column with general C, Chl, and nutrient concentrations (see
Table A3.3). Advection and diffusion rates, as well as temperature, for each vertical grid box were
averaged through time over the entire simulation period for each trajectory and the mean was used
for the spin-up period. For each latitude where particles were released, the biogeochemical initial
conditions were determined using physical values from the most offshore particle location at that
latitude. The biogeochemical model was run for 21 days to allow all variables to reach steady
state. The steady-state values were then saved and used as the initial conditions for all particles
released at the corresponding latitude. Sensitivity analyses show that the impact of the chosen
initial conditions were negligible after the first 10 days of the simulation (Figure A3.2). Therefore,
in the following analyses, we omit the first 10 days to discard any potential impacts from the chosen
initial conditions.
4.2.5 Post-processing model analyses
In subsequent analyses comparing the different rates of limitation (temperature, nutrient, and light),
we used a normalized and filtered version of the data to remove the impact of the diurnal cycle.
A 24-hour moving mean was applied to the data to remove the impact of night time light limitation. Additionally, at each time step, the light limitation values calculated in Equation 13 were
normalized through the whole water column by the maximum value (least limited) within the water
75
column at that time step. This value was typically located at the surface but due to photoinhibition
in the model, could sometimes be slightly below the surface level.
4.3 Results
Our framework allowed us to quantify bio-physical interactions occurring along a Lagrangian trajectory. Critically, we resolved the critical temporal and spatial scales to accurately model the
biological response to fine-scale physics. Below we present the validation of the model dynamics,
an assessment of the controls on biomass accumulation in the Lagrangian reference frame, and the
impact of environmental variability on biomass accumulation.
4.3.1 Model Validation
The results from the biogeochemical model were validated using the extensive California Cooperative Oceanic Fisheries Investigations (CalCOFI) observational dataset, as well as measurements
from Li et al., 2010. The CalCOFI program has been sampling the study region, which spans
from Southern California to San Francisco, quarterly since 1949. We compared the model output
(March 27 - June 1) to the all measurements performed in the month of April from CalCOFI. Li et
al., 2010 sampled the upwelling period (May) in 2006 which corresponds well with the seasonality
simulated by the model. Model data used for validation were extracted from our water column
model at the depths occupied by the Lagrangian drifters at each time step.
4.3.1.1 Nitrate Concentration
The Southern California Current system is classically a nutrient limited system with low concentrations of nutrient in the surface waters (Deutsch et al., 2021; Guo et al., 2014; Guo et al., 2020). This
nutrient limitation is relieved in the spring with upwelling of deep waters which bring nutrients up
into the sunlit surface. The model captures the observed profile of NO3 with depth (Figure 3).
Surface (<10m) nitrate in the model was on average 2.5 mmol/m3 with a range of 0-6.4 mmol/m3
(2σ) which is consistent with the CalCOFI data (average NO3 of 4.4. mmol/m3 with a range of
76
0-17.1 mmol/m3
). Li et al. (2010) had significantly fewer samples (N=38 vs N=1.9 million for the
model and N=10,183 for CalCOFI bottle data) but showed a similar range for surface NO3. The
increase in NO3 with depth in the model was also consistent with both the CalCOFI and the Li et
al. (2010) data (Figure 3). Specifically, NO3 increased with depth in the model at an average rate
of is -0.12 ±0.004 (2σ) m/mmol NO3/m3
consistent with the CalCOFI nitrate data which increased
with depth at an average rate -0.11 ±0.009 (2σ) m/mmol NO3/m3
. The measurements from Li et
al. (2010) also fit within the overall range of the model and bottle data. These results instill confidence that the physical and biological processes that control the vertical distribution of nitrate are
well represented in our model.
Figure 4.3: The relationship between nitrate and depth for the CalCOFI bottle data (light blue), our model
data (dark blue), and measurements from Li et al. (2010) (orange diamonds). The yellow and blue lines represent the bottle and model means, respectively. Dashed lines of the same colors represent the 2σ standard
deviation.
In the Southern California region of the CCS, nitrate concentrations largely follow a linear
relationship with temperature within 12 km of the shore (Lucas et al., 2011b). The relationship
between nitrate and temperature in the CalCOFI dataset, the Li et al. (2010) measurements, and our
77
water column model all followed a similar linear relationship (Figure 4). Generally, as temperature
decreased, nitrate concentration increased. Overall, the rate of draw down of nitrate as waters
warm was similar to that observed in the CalCOFI data in waters <12°C with the rate of nitrate
decrease in the CalCOFI data of -6.04 mmol NO3/m3
/◦C and a rate in the model of -1.48 mmol
NO3/m3
/◦C. However, the model underestimated NO3 concentrations in cooler waters and over
estimated concentrations in warmer waters (Figure 4). Given the consistency in the depth profiles
of nitrate between the model and the CalCOFI data both in terms of the mean but also in the spread
of the NO3 data, we hypothesize that the mismatch in the nitrate versus temperature relationship
may in part be due to under-sampling by the model of the cold (deep), nitrate rich waters and
over-sampling of the lower nitrate, surface waters.
Figure 4.4: The relationship between nitrate and temperature for the bottle date, the Li et al., 2010 measurements, our model data, and the Lucas model (Lucas et al., 2011). Solid dark blue points and light blue
points are the model and CalCOFI bottle data, respectively. The yellow and blue lines represent the bottle
and model means, respectively. Dashed lines of the same colors represent the 2σ standard deviation. Orange
diamonds represent the measurements from Li et al., 2010. The gray solid line shows the temperature-nitrate
relationship Lucas et al., 2011b.
78
4.3.1.2 Light, Photosynthesis, and Growth
Phytoplankton growth is limited by nutrients, temperature, and light. As we demonstrated above,
the model captured the observed depth profile of nitrate. We also found good agreement with light
attenuation in our model compared with measurements of photosynthetically available radiation
(PAR) taken in the CCS (Mitchell et al., 2017) (Figure 5). The emergent light field in the model
shows a strong and significant correlation (R2=1, p<0.05) with the mean of the measured data and
the 2σ confidence interval encompasses nearly all the data points from Mitchell et al., 2017.
0 20 40 60 80 100
% Surface PAR
-110
-100
-90
-80
-70
-60
-50
-40
-30
-20
-10
0
Depth [m]
Mitchell et al., 2017
model mean
Figure 4.5: Vertical profile of mean PAR in the model compared to observed PAR data from Mitchell et al.,
2017.
This investigation examines the influence of fine-scale bio-physical interactions on phytoplankton growth and biomass accumulation. Chlorophyll, the crucial link between light availability and
growth rates, exhibits a distinct vertical distribution. Observational data reveal relatively elevated
chlorophyll concentrations near the surface, with a progressive decline with increasing depth (Figure 6), a pattern that is successfully reproduced by the model. Overall, the model represents the
79
two types of responses observed in the Li et al. data - elevated surface Chl with a maximum around
6 mg Chl/m3
and a more uniform distribution of Chl with depth. The model Chl concentrations
also generally agree with the observed Chl pattern in the CalCOFI bottle data but the model consistently overestimates Chl by approximately 1 mg Chl/m3
. The model mean closely tracks estimated
growth rates from Li et al., (2010) (Figure 7). Additionally, 95% of the Li et al., (2010) growth
rates were within the 2σ standard deviation of the model mean. Similarly, the model captured the
observed vertical profile for carbon reported by Li et al. (Figure 8). The overall good agreement
between the model and observational data for carbon, chlorophyll, and growth rates suggests that
our model captures the primary processes driving these dynamics in the CCS.
Figure 4.6: Vertical profiles of Chl with depth. Data from CalCOFI (light blue dots, yellow lines), our
model (dark blue dots, green lines), and Li et al. (2010) (orange diamonds) are shown.
80
Figure 4.7: Realized growth rates from the model (dark blue dots) compared to estimated growth rates from
Li et al. (2010) (orange dots) as a function of depth. The model mean( green line) and standard deviation
(dotted lines, 2σ) encompass 95% of the Li et al. (2010) data.
81
Figure 4.8: Vertical profiles of carbon from the model (dark blue dots) and Li et al. (2010). The model
captures both the extreme and general profile of carbon from the Li et al. data.
82
4.3.2 Biological and Physical Control on Biomass Accumulation
Coupling Lagrangian particle trajectories with the biogeochemical model allowed us to understand
biological and physical drivers of phytoplankton biomass accumulation in the CCS. In general,
tight coupling between growth and loss limited the accumulation of biomass. Thus, it was only
where these processes become decoupled that we observed significant increases in phytoplankton
biomass. Here, we deconvolved the factors controlling changes in biomass to identify the relative
contributions of biological (growth, mortality) and physical (dilution from diffusive mixing) processes. Specifically, for each Lagrangian trajectory, we identified the primary processes driving
changes in phytoplankton biomass at each time step as the biological or physical processes that
drove the accumulation or loss in biomass. The changes in biomass due to the diel cycle were
filtered out using a 24-hour moving mean for all variables.
Overall, we observed that physical processes and biological processes had approximately equal
impact on changes in biomass in the CCS, with biology dominating the change in biomass on
average 44% of the time (Figure 9). As expected, biological processes resulted in biomass accumulation 99% of the time as a result of phytoplankton growth (Figure 10). Biomass loss due to
biological processes occurred below the euphotic zone when light limitation caused growth rates
to be zero and respiration and mortality resulted in overall biomass loss. Conversely, when physical processes dominated biomass change a loss in biomass was observed the majority (84%) of
the time (Figure 10). This loss primarily occurred due to diffusive mixing which acted to dilute
biomass concentrations.
4.3.2.1 Biological Controls
To grow, phytoplankton must overcome chemical (i.e., nutrient) and physical (i.e., temperature
and light) limitations on growth (see Eq 12-17). Along the Lagrangian trajectories, when biomass
was accumulating, nutrients were the least limiting factor the majority of the time (73% of the
time, Figure 11a). This suggests that most of the biomass accumulation in the CCS occurs when
nutrient limitation is removed, as opposed to alleviating light or temperature limitation. This is
83
0 200 400 600 800 1000 1200 1400
Lagrangian Trajectory Number
0
0.2
0.4
0.6
0.8
1
Proportion of time controlling biomass
light lim
nutrient lim
temp. lim
diffusion
advection
Figure 4.9: Proportion of time that changes in biomass are controlled by biological and physical processes.
Data are ordered from highest proportion of biology dominated changes to lowest. Biological controls
on biomass change are light limitation, nutrient limitation, and temperature limitation. When biological
processes control the change in biomass, the majority of that time nutrient limitation is the dominant factor.
Physical processes that impact biomass concentration are advection and diffusion. Diffusion is the most
dominant physical process.
consistent with our expectation that this is primarily a nutrient limited system. Additionally, when
growth could not outpace mortality (biomass was decreasing), light limitation and nutrients impacted accumulation nearly equally, though light limitation was slightly more common (54% of
the time, Figure 11b). This suggests that biologically mediated biomass decreases occurred either
when trajectories were both subducted out of the euphotic zone or when nitrate uptake outpaced
replenishment in the euphotic zone. Indeed, the average depth for time-steps where biomass was
decreasing and light was most limiting was 56.3 m while when biomass was decreasing but nutrients were more limiting than light, the average particle depth was 8.4 m.
84
-0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6
ΔBiomass [mmol C/m³/hr]
0
1
2
3
4
5
6
Number of timesteps
10 4
physics
biology
phys mean
bio mean
Figure 4.10: When biology was the controlling process, the change in biomass was on average positive
(0.18 mmol C/m3
/hr or 4.32 mmol C/m3
/day). When diffusive mixing was the controlling process, the
change in biomass skewed negative with an average change of -0.19 mmol C/m3
/hr (4.6 mmol C/m3
/day)
85
Most limiting factor when ∆biomass < 0
0 100 200 300 400 500
Lagrangian Trajectory Number
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Proportion of growth rate limitation
light
temperature
nutrients
Least limiting factor when Δbiomass > 0
0 200 400 600 800 1000 1200 1400
Lagrangian Trajectory Number
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Proportion of growth rate limitation
light
temperature
nutrients
a b
Figure 4.11: Limiting biological factors when changes in biomass are due to biological processes. (a)
When biomass concentration increases due to biological processes, nutrients are typically the least limiting
factor, although relieving light limitation is also important. (b) When biomass concentration decreases due
to biological processes, it is approximately half the time due to light limitation while the other half is due to
nutrient limitation.
86
4.3.2.2 Physical Controls
One of the key advantages of our Lagrangian-based model was the ability to track the impact of
diffusive mixing on biomass accumulation. Diffusive mixing played a key role in biomass accumulation in the majority of our trajectories. Specifically, for all particle timesteps in which the
biomass concentration was primarily changed via physical processes, the change was overwhelmingly due to diffusive mixing (Figure 12). Nearly all (99%) of the timesteps in which biomass
decreased was due to mixing of high biomass surface waters with low biomass deep waters (dilution effect). The mean particle depth at which diffusive mixing decreased biomass was 15.0m.
This mixing also acted to decrease the mean light level available for growth thereby further reducing biomass accumulation. Similarly, turbulent mixing was the primary process (80% of the
time) when physical dynamics resulted in increased biomass. This occurred when the particle was
located deeper in the water column but above the euphotic zone depth. Typically in this scenario
the base of the mixed layer was below the base of the euphotic zone and mixing brought biomass
from the surface downward. The average particle depth that experienced an increase in biomass
due to diffusive mixing was 33.4m. Given the Lagrangian framework of our model, we expected
advective losses to be minimal as we were evaluating changes along an advected path. Consistent
with this, biomass changes due to advection accounted for <20% of all biomass changes due to
physics (Figures 12 and 9).
4.3.2.3 Critical Diffusivity
The critical diffusivity (Kc) is defined as the mixing rate at which phytoplankton growth can overcome the dilution effect induced by mixing and biomass can accumulate (Taylor et al., 2012). Previous work proposed that in the North Atlantic, the spring phytoplankton bloom (positive biomass
accumulation) occurs as the turbulent mixing characteristic of winter gives way to a calmer, more
stable, spring time water column (Mahadevan et al., 2012; Taylor et al., 2012; Taylor, 2016). The
North Atlantic is primarily a light limited system (e.g. Lacour et al., 2015; Rumyantseva et al.,
2019) unlike the CCS. However, some previous work in the CCS suggested that water column
87
Physical controls when Δbiomass > 0
0 200 400 600 800 1000
Lagrangian Trajectory Number
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Proportion of total change due to physics
diusion
advection
Physical controls when Δbiomass < 0
0 200 400 600 800 1000 1200
Lagrangian Trajectory Number
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Proportion of total change due to physics
diusion
advection
Figure 4.12: Dominating physical process when the change in biomass is a) positive, b) negative. Turbulent/diffusive mixing is the controls biomass most often but advection does sometimes lead to biomass
increases.
stabilization for at least 3 days, potentially up to 7 days, was necessary in order for biomass to accumulate (Dugdale et al., 1989; Wilkerson et al., 2006). The Lagrangian framework of our model
allowed us to directly track the combinations of mixing and growth rates and relative timescales
that permit biomass accumulation. We compared the change in biomass per time step to the diffusivity and growth rate of the particle. Additionally, we were able to record the time in between
mixing events along the trajectories and compare these timescales with the modeled changes in
biomass.
The observed relationship between biomass accumulation and diffusive mixing highlights the
nuanced complexity of bio-physical interactions in the CCS. Rather than a simple relationship of
growth rate outpacing dilution, the model provides insight into the dynamic and complicated relationships between mixing, nutrient limitation, light limitation, and biomass accumulation rate.
Consistent with the critical diffusivity theory, the highest biomass accumulation rates were associated with the lowest diffusive mixing rates (Figure 13). However, low mixing rates (Kz <2e-3
m2
/s) were sometimes associated with biomass losses (30% of the time) even when growth rates
were significant. As diffusive mixing increased, we observed an average decrease in biomass accumulation rates with a little over half (≈54%) of time-points resulting in biomass loss when mixing
rates were between 5e-3 m2
/s and 15e-3 m2
/s. Unexpectedly, the time steps within this range also
88
Figure 4.13: Relationship between diffusivity (Kz) and the change in biomass per time step. Colors of the
data points represent the growth rate. Particles with the highest growth rates are counter intuitively most
associated with biomass decreases.
had the highest growth rates. As mixing increased above 20e-3 m2
/s, there was no clear relationship between mixing and biomass accumulation. Specifically, in direct contrast with the critical
diffusivity hypothesis, these high mixing regimes were nearly equally associated with biomass accumulation (42% of timesteps) and biomass loss (58%). Under higher mixing regimes (>20e-3
m2
/s), growth rates were substantially lower 0.42 day-1 ± 0.43 day-1 (2σ) than compared to 0.54
day-1 ± 0.63 day-1 when mixing was between 5e-3 m2
/s and 15e-3 m2
/s. Overall, our analysis
suggests that for the CCS mixing is not the dominant process impacting biomass accumulation
at these high mixing levels – at least when evaluating bio-physical interactions in the Lagrangian
framework.
4.3.3 The Impact of Lagrangian Variability
To demonstrate the effect of Lagrangian variability on biomass accumulation, we considered a set
of Lagrangian particles that were initialized 395 km apart but ended the simulation 3.1km away
from each other and at similar depths (Figure 14, 13.8m for Example Trajectory 1 and 6.0m for
Example Trajectory 2). These two trajectories provide an illustrative example of the importance
89
of Lagrangian variability. Specifically, both trajectories had similar final temperatures – 14.4 for
Example Trajectory 1 (ET1) and 14.8 for Example Trajectory 2 (ET2) – and similar mean environmental conditions over the trajectory: temperature (11.8 - ET1, 11.5 - ET2), nitrate concentration
(5.9 mmol/m3
- ET1, 5.1 mmol/m3
- ET2), mixed layer depth (31.8m - ET1, 41.0m - ET2), and
photic zone depths (66.8m - ET1, 64.5m - ET2). Additionally, the particles encountered similar
patterns in diffusive mixing with similar timing of increases and decreases in mixing rates (Figures
15).
-130 -129 -128 -127 -126 -125 -124
Longitude
35.5
36
36.5
37
37.5
38
38.5
39
39.5
40
40.5
Latitude
Example Trajectory 1
Example Trajectory 2
Figure 4.14: Two example trajectories over the 66.7 day simulation. Stars mark the start of the simulation.
These particles started 395km apart and ended only 3.1km from one another.
Despite ending in the same geographical location and having largely similar mean environmental conditions, ET1 has 54% more biomass than ET2 in the final timestep (Figures 15a, 15b).
This difference was primarily driven by differential physics experienced by the two drifters over
90
the last 23 days of the simulation. Around day 43, ET1 was swiftly subducted 38.3m over 1.4
days and remained below the mixed layer and above the euphotic zone for 11 days. Because the
particle along ET1 was within the euphotic zone, the growth rate was positive but low (0.1 day−1
on average) due to severe light limitation (γL = 0.02 on average) During this time, strong downward velocities advected biomass from the surface waters, resulting in an increase in subsurface
biomass at the depth of the ET1 particle.
Conversely, during this period, ET2 remained at the surface within both the mixed layer and
euphotic zone . The mean growth rate for ET2 over this period was substantially higher, 0.7 day−1
,
as the particle was not severely limited by nutrients, light, or temperature. However, variable
physics in particular periodic elevated diffusive mixing, acted to dilute the biomass along ET2 thus
preventing substantial biomass accumulation during this period. The net result of these processes
was that by day 66.7, the particle along ET1 had a biomass of 23.7 mmol C/m3 whereas ET2 had
a final biomass concentration of 15.4 mmol C/m3
.
-200
-180
-160
-140
-120
-100
-80
-60
-40
-20
0
Depth [m]
0
5
10
15
20
25
Biomas [mmol C m-3]
Biomass
Particle
MLD
Photic Zone
0
5
10
15
20
25
Biomass [mmol C m-3]
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.1
Kz [m2⁄s]
Carbon
Mixing Rate
EXAMPLE TRAJECTORY 1
10 20 30 40 50 60
Time [days]
-200
-180
-160
-140
-120
-100
-80
-60
-40
-20
0
Depth [m]
0
5
10
15
20
25
30
35
Nitrate
Particle
MLD
Photic Zone
0
5
10
15
20
25
Biomass [mmol C m-3]
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.1
Kz [m2⁄s]
Carbon
Mixing Rate
EXAMPLE TRAJECTORY 2
10 20 30 40 50 60
Time [days]
-200
-180
-160
-140
-120
-100
-80
-60
-40
-20
0
Depth [m]
0
5
10
15
20
25
30
35
NO3 [mmol⁄m3}
Nitrate
Particle
MLD
Photic Zone
NO3 [mmol⁄m3}
a b
f
c
e
10 20 30 40 50 60
time [days]
-200
-150
-100
-50
0
Depth [m]
0
5
10
15
20
25
Biomass [mmol C/m3]
Biomass
Particle
MLD
Photic Zone
d
Figure 4.15: Comparison of two example trajectories. (a,b) Biomass concentration through time on the left
axis and on the right axis are the diffusive mixing rates. Shaded regions represent timesteps in which the
change in biomass was dominated by physical processes. (c,d) The biomass concentration along the entire
200m water column for the whole simulation period. Also shown are the depths of the particle along the
trajectory, the mixed layer, and the base of the euphotic zone. (e,f) Same as in (c,d) but nitrate concentration
instead of biomass.
91
0 5 10 15 20 25
ET1 biomass [mmol C/m3]
0
5
10
15
20
25
ET2 Biomass [mmol C/m3]
-200
-180
-160
-140
-120
-100
-80
-60
-40
-20
0
Depth [m]
0 5 10 15 20 25 30
ET1 Nitrate [mmol/m3]
0
5
10
15
20
25
30
ET2 Nitrate [mmol/m3]
-200
-180
-160
-140
-120
-100
-80
-60
-40
-20
0
Depth [m]
a b
Figure 4.16: Comparison of the whole water column nitrate (a) and biomass (b) concentrations of ET1 and
ET2 at the final timestep of the simulation.
To further support that the differences in biomass were due to Lagrangian history and not that
the trajectories ended at slightly different depths, we compared the biomass and nitrate concentrations at the final time step for the whole water column model (200m). ET1 consistently had higher
biomass and lower nitrate than ET2 (Figure 16) by up to -31% for nitrate and 41% for biomass. The
top and bottom of the two water column models were the same as those were set by the boundary
conditions. The differences in biomass between the two particles reflect differences in the entire
water column.
4.4 Discussion
Our novel modeling approach allowed us to, for the first time, investigate how fine-scale physical
dynamics impact phytoplankton growth along Lagrangian trajectories. We used a high-resolution
(500 m), physical oceanography model of the CCS and Lagrangian particle tracking paired with
a biogeochemical model to assess the impact of Lagrangian biophysical interactions on biomass
accumulation. We found that Lagrangian variability impacts biomass accumulation in two ways.
First, in our data, a critical diffusivity value, over which biomass cannot accumulate, did not exist. All mixing rates were associated with both biomass accumulation and loss. Second, our data
92
demonstrated that within the Lagrangian reference frame there was considerable variability between trajectories. Although two particles end at the same location, they do not have the same
biomass concentration due to experiencing different environmental histories.
Previous modeling work in the CCS suggested that rapid subduction due to fine-scale physics,
particularly submesoscale eddies and fronts, quickly removes phytoplankton (and therefore biomass)
out of the photic zone and limits growth rates and nitrate uptake relative to coarse grid models
(Kessouri et al., 2020; Damien et al., 2023). Our results indicate that while biomass loss due to
light limitation does happen, it is rare and occurs only 1% of the time that biomass change is dominated by biology. Rather, our results suggest that nearly all biomass losses are due to diffusive
mixing. Additionally, we have shown that while Lagrangian particles in the CCS do experience
rapid vertical changes (Chapter 2, see also Figure 15 here), we also found that Lagrangian particles do not encounter submesoscale features frequently nor interact with the features for very long.
Rather, we propose that physical dynamics that impact mixing would play a larger role than submesoscale features. One such possibility are internal gravity waves. When internal gravity waves
break, they also cause mixing that can increase nutrient concentrations through mixing between
surface and deep waters (Sandstrom et al., 1984; Lucas et al., 2011a) or potentially dilute surface
biomass. It has been previously suggested that internal gravity waves may cause phytoplankton
patchiness (Kamykowski, 1974) through modulating light limitation via vertical motion (Evans
et al., 2008; Garwood et al., 2020). For the CCS, the vertical motions are dominated by internal
gravity waves which are associated with similar vertical velocities and timescales as submesoscale
eddies and fronts (Torres et al., 2022). Our work expands on this previous study and suggests that
internal gravity waves are the dominant fine-scale processes impacting phytoplankton growth and
biomass accumulation in the CCS.
We also hypothesize that better resolved diffusive mixing would result in a decrease in biomass
loss relative to a coarse-grid model, as was found in Kessouri et al., 2020; Damien et al., 2023.
As model resolution increases, smaller-scale processes (Trotta et al., 2017), specifically mixing
(Liu et al., 2021) start to emerge. Our model results show that even small increases in mixing
93
rates can result in biomass loss (Figure 15). By increasing the horizontal resolution of physical
ocean models, diffusive mixing becomes more apparent. When comparing two models of different
horizontal resolutions, decreases in surface biomass could be due, in part, to an increase in dilution
via diffusive mixing.
In addition to the physical processes, we show that overcoming nitrate limitation is also required for biomass accumulation. However, one limitation of our model is that is does not account
for iron, a micronutrient that can be limiting in parts of the CCS (Kudela et al., 2008). Incorporating iron dynamics into future models could provide a more comprehensive understanding of the
factors controlling phytoplankton growth in this region.
4.5 Conclusions
Our high-resolution model and Lagrangian drifter tracking approach provided insight into the complex biophysical interactions governing phytoplankton biomass accumulation in the CCS. We highlight the role of turbulent mixing in biomass accumulations rates, challenging previous assumptions about the importance of rapid subduction. Establishing the influence of turbulent mixing
versus vertical advection would not have been possible to determine without modeling in the Lagrangian reference frame. With our novel framework, we were able to determine at each timestep
the dominant controls on phytoplankton growth and biomass accumulation.
Future research should aim to bridge the gap between Lagrangian and Eulerian frameworks
to better interpret the impact of environmental variability. More Lagrangian-based studies that
incorporate interannual variations in upwelling will also provide better insight into complimentary roles of submesoscale features and turbulent mixing. By addressing these knowledge gaps,
we can further our understanding of the relationship between physical processes and phytoplankton productivity in the CCS, contributing to more accurate predictions of ecosystem responses to
environmental change.
94
Chapter 5
Conclusions and Future Work
5.1 Conclusions
Phytoplankton are microorganisms that have a global impact on biogeochemical cycles. With improved horizontal grid resolution, the inclusion of phytoplankton dynamics in large scale biogeochemical models is also improving. However, due to computational constraints, global simulations
do not resolve environmental variability occurring at the sub-grid scales that is crucial to phytoplankton dynamics. Additionally, simplifying assumptions about the succession of phytoplankton
phenotypes under changing environmental conditions are also required to lower computational
requirements which may not reflect actual phytoplankton community dynamics. By understanding fine-scale phytoplankton growth dynamics better, we can begin to appropriately account for
small-scale processes in large-scale models.
In this body of work, I have demonstrated that environmental variability differs between the
Lagrangian and Eulerian reference frames. Specifically, a passive drifter will encounter more
temperature variability than a fixed point in the Southern Ocean. However, the differences in
environmental variability between reference frames found in the Southern Ocean cannot be generalized to all locations. In the California Current System, the difference between Eulerian and
Lagrangian environmental variability depended on direction of Lagrangian particle advection and
the variable of interest (e.g. temperature or mixing rates). Specifically, as a particle moves offshore
it encounters more temperature variability than the corresponding Eulerian reference frame while
the opposite is true when a particle moves onshore. Conversely, as a particle move offshore it
encounters less mixing variability than the corresponding Eulerian reference frame and vice versa
for particles advected onshore.
95
Standard methods for modeling phytoplankton growth dynamics do not accurately capture the
response of biomass accumulation to environmental variability because the standard parameterization does not account for the timescale of response from the phytoplankton community. One
of the simplifying assumptions made within most models is that within a phytoplankton community, there will be a smooth transition of the phenotype best suited in the current environment to
the best suited in the next environment. The problem with this is twofold. First, there is a lag
between when the temperature changes and when the biomass of the best suited individual accumulates sufficiently to impact overall community growth rates. This lag is a function of the rate
and magnitude with which temperature changes. Second, due to the variable nature of the natural environment, conditions are constantly changing, typically faster than the rate of community
adjustment, and therefore the phytoplankton community composition is rarely in equilibrium with
the instantaneous environmental conditions. The time required for a community to reach equilibrium will depend on the rate of acclimation of each individual phenotype within the community
and will combine non-linearly with variability in other environmental conditions such as light and
nutrient concentration.
With each iteration of increasing the horizontal resolution in biogeochemical models, the impact of smaller-scale physical dynamics on phytoplankton biomass accumulation becomes clearer.
The results presented in this dissertation demonstrate that changes in phytoplankton biomass that
occur when submesoscale dynamics are resolved are not specifically related to the eddies and
fronts that emerge. Rather, changes in the mixed layer dynamics as a result of increased horizontal
resolution had the largest impact on biomass concentration. As model resolution continues to increase and smaller-scale dynamics are fully resolved, large-scale dynamics are also better resolved.
The results in this dissertation suggest that the increased ability to resolve mixed-layer dynamics
because of better resolved energetic feedback between small- and large-scale dynamics will have
the most impact on biomass concentration, not specifically the smaller mechanisms themselves.
While global models currently lack the computational power to fully resolve submesoscale dynamics, their impact on mixed layer processes can be parameterized. High-resolution regional models
96
can help establish relationships between mixed layer depth variability and phytoplankton biomass,
enabling the incorporation of these parameterizations into coarser global models. This approach
offers a way to account for the influence of unresolved mixed layer dynamics on phytoplankton,
even in computationally limited global simulations.
Overall, I demonstrated that spatial and temporal heterogeneity largely impact overall phytoplankton biomass accumulation dynamics in a way that is not captured by most modeling techniques. However, I also showed the potential for including the impact of environmental variability
on biomass accumulation in biogeochemical models through establishing relationships between
what phytoplankton experience and the conditions simulated.
5.2 Future Work
Future work should strive to identify overarching principles that govern the interplay between
Lagrangian and Eulerian environmental variability. Establishing these connections, particularly
for biologically important variables like temperature, light, nutrients, and mixing, is essential for
successfully incorporating the impact of Lagrangian processes into the framework of global models. In addition to considering Lagrangian environmental variability, we must also factor in the
timescales of phytoplankton response. While Eulerian models can reproduce observed spatial patterns in phytoplankton biology, this alone is insufficient for a comprehensive understanding of the
bio-physical dynamics governing global biogeochemical cycling. This limitation arises because of
the lag-time in biomass response to environmental changes. The biomass accumulation at a fixed
Eulerian location reflects not just local conditions but also the influence of upstream dynamics.
For instance, spatial patterns in phytoplankton biology can arise from either phytoplankton
community growth or surface convergence, each with distinct implications for nutrient cycling
and photosynthesis rates. Without a Lagrangian perspective that tracks the movement and history
of water parcels, distinguishing between these two processes becomes challenging, hindering our
ability to fully grasp the complex interplay of factors shaping marine ecosystems.
97
Future observational work would benefit from higher resolution sampling in space, and time,
but also from sampling taken in a truly Lagrangian reference frame. Autonomous vehicles, e.g.
wire-walkers, are limited to the influence of only surface currents and will inherently over-sample
regions of surface convergence. Biogeochemical floats provide a great method to sample the water
column as a whole but are only quasi-Lagrangian due to the parking depths and vertical transects they take. Ideally, we would develop a device that can be set to be neutrally buoyant at
specified densities that would then advect vertically and horizontally. By deploying a network of
these types of floats in a three-dimensional release pattern, not only would they track Lagrangian
environmental variability, but they could have sensors onboard that would track important biogeochemical parameters such as light level, turbidity, temperature, salinity, and particle size. Further,
by employing a network of these particles, the relative motions of the particles could be used
to calculate divergence which could potentially aid in submesoscale feature identification along
these Lagrangian trajectories. From these hypothetical, Lagrangian datasets, we could then get a
more realistic representation of environmental variability and incorporate that into biogeochemical
models.
In conclusion, the findings of this dissertation reveal that the intricate interplay between spatial
and temporal heterogeneity significantly influences phytoplankton biomass accumulation, in ways
that are not adequately captured by current modeling techniques. By establishing discrepancies
between Lagrangian and Eulerian reference frames, challenging traditional modeling assumptions,
and delving into the effects of submesoscale dynamics, this work underscores the limitations of
existing biogeochemical models and paves the way for future refinements. While phytoplankton
responses to environmental variability are undeniably complex, the potential to enhance biogeochemical models through improved parameterizations and the integration of Lagrangian perspectives is evident. By bridging the gap between the lived experience of phytoplankton and the simulated world of models, we can strive towards a more accurate and holistic understanding of global
biogeochemical cycling - one that not only reproduces spatial patterns but also illuminates the
underlying bio-physical mechanisms driving them.
98
References
Abbott, Mark R. et al. (1990). “Observations of Phytoplankton and Nutrients from a Lagrangian
Drifter off Northern California”. In: Journal of Geophysical Research: Oceans 95.C6, pp. 9393–
9409. ISSN: 2156-2202. DOI: 10.1029/JC095iC06p09393.
Artana, Camila et al. (2019). “The Malvinas Current at the Confluence With the Brazil Current: Inferences From 25 Years of Mercator Ocean Reanalysis”. In: Journal of Geophysical Research:
Oceans 124.10, pp. 7178–7200. ISSN: 2169-9275, 2169-9291. DOI: 10.1029/2019JC015289.
Arteaga, Lionel A. et al. (2020). “Seasonal Modulation of Phytoplankton Biomass in the Southern
Ocean”. In: Nature Communications 11.1, p. 5364. ISSN: 2041-1723. DOI: 10.1038/s41467-
020-19157-2.
Balwada, Dhruv et al. (2021). “Vertical Fluxes Conditioned on Vorticity and Strain Reveal Submesoscale Ventilation”. In: Journal of Physical Oceanography 51.9, pp. 2883–2901. ISSN:
0022-3670, 1520-0485. DOI: 10.1175/JPO-D-21-0016.1.
Barton, Andrew D. et al. (2010). “Patterns of Diversity in Marine Phytoplankton”. In: Science
327.5972, pp. 1509–1511. ISSN: 0036-8075, 1095-9203. DOI: 10.1126/science.1184961.
Baudry, Jer´ emy et al. (2018). “Turbulent Mixing and Phytoplankton Life History: A Lagrangian ´
versus Eulerian Model Comparison”. In: Marine Ecology Progress Series 600, pp. 55–70. ISSN:
0171-8630, 1616-1599. DOI: 10.3354/meps12634.
Beal, Lisa M. et al. (2015). “Capturing the Transport Variability of a Western Boundary Jet: Results
from the Agulhas Current Time-Series Experiment (ACT)”. In: Journal of Physical Oceanography 45.5, pp. 1302–1324. ISSN: 0022-3670. DOI: 10.1175/JPO-D-14-0119.1.
Bernhardt, Joey R et al. (2018). “Nonlinear Averaging of Thermal Experience Predicts Population
Growth Rates in a Thermally Variable Environment”. In: Proceedings of the Royal Society B:
Biological Sciences 285.20181076, p. 10.
Bopp, L. et al. (2013). “Multiple Stressors of Ocean Ecosystems in the 21st Century: Projections
with CMIP5 Models”. In: Biogeosciences 10.10, pp. 6225–6245. ISSN: 1726-4189. DOI: 10.
5194/bg-10-6225-2013.
Boyd, Philip W. (2019). “Physiology and Iron Modulate Diverse Responses of Diatoms to a Warming Southern Ocean”. In: Nature Climate Change 9.2, pp. 148–152. ISSN: 1758-678X, 1758-
6798. DOI: 10.1038/s41558-018-0389-1.
Boyd, Philip W. et al. (2013). “Marine Phytoplankton Temperature versus Growth Responses from
Polar to Tropical Waters – Outcome of a Scientific Community-Wide Study”. In: PLoS ONE
8.5. Ed. by Howard Browman, e63091. ISSN: 1932-6203. DOI: 10 . 1371 / journal . pone .
0063091.
Boyd, Philip W. et al. (2016). “Biological Responses to Environmental Heterogeneity under Future
Ocean Conditions”. In: Global Change Biology 22.8, pp. 2633–2650. ISSN: 1365-2486. DOI:
10.1111/gcb.13287.
Brannigan, L. (2016). “Intense Submesoscale Upwelling in Anticyclonic Eddies”. In: Geophysical
Research Letters 43.7, pp. 3360–3369. ISSN: 1944-8007. DOI: 10.1002/2016GL067926.
Brock, Thomas D. (1981). “Calculating Solar Radiation for Ecological Studies”. In: Ecological
Modelling 14.1, pp. 1–19. ISSN: 0304-3800. DOI: 10.1016/0304-3800(81)90011-9.
99
Buitenhuis, Erik T. et al. (2013). “Combined Constraints on Global Ocean Primary Production
Using Observations and Models”. In: Global Biogeochemical Cycles 27.3, pp. 847–858. ISSN:
1944-9224. DOI: 10.1002/gbc.20074.
Busecke, Julius J. M. et al. (2019). “The Equatorial Undercurrent and the Oxygen Minimum Zone
in the Pacific”. In: Geophysical Research Letters 46.12, pp. 6716–6725. ISSN: 1944-8007. DOI:
10.1029/2019GL082692.
Caldeira, R. M. A. et al. (2005). “Island Wakes in the Southern California Bight”. In: Journal of
Geophysical Research: Oceans 110.C11. ISSN: 2156-2202. DOI: 10.1029/2004JC002675.
Calil, Paulo H. R. et al. (2021). “Filaments, Fronts and Eddies in the Cabo Frio Coastal Upwelling
System, Brazil”. In: Fluids 6.2, p. 54. ISSN: 2311-5521. DOI: 10.3390/fluids6020054.
Cao, Haijin et al. (2024). “Isopycnal Submesoscale Stirring Crucially Sustaining Subsurface Chlorophyll Maximum in Ocean Cyclonic Eddies”. In: Geophysical Research Letters 51.4, e2023GL105793.
ISSN: 1944-8007. DOI: 10.1029/2023GL105793.
Capet, X. et al. (2008a). “Mesoscale to Submesoscale Transition in the California Current System.
Part I: Flow Structure, Eddy Flux, and Observational Tests”. In: Journal of Physical Oceanography 38.1, pp. 29–43. ISSN: 0022-3670, 1520-0485. DOI: 10.1175/2007JPO3671.1.
— (2008b). “Mesoscale to Submesoscale Transition in the California Current System. Part II:
Frontal Processes”. In: Journal of Physical Oceanography 38.1, pp. 44–64. ISSN: 0022-3670,
1520-0485. DOI: 10.1175/2007JPO3672.1.
Castelao, Renato M. et al. (2006). “Sea Surface Temperature Fronts in the California Current System from Geostationary Satellite Observations”. In: Journal of Geophysical Research: Oceans
111.C9. ISSN: 2156-2202. DOI: 10.1029/2006JC003541.
Centurioni, L. R. et al. (2008). “Permanent Meanders in the California Current System”. In: Journal of Physical Oceanography 38.8, pp. 1690–1710. ISSN: 0022-3670, 1520-0485. DOI: 10.
1175/2008JPO3746.1.
Chavez, Francisco P. et al. (1991). “Horizontal Transport and the Distribution of Nutrients in the
Coastal Transition Zone off Northern California: Effects on Primary Production, Phytoplankton Biomass and Species Composition”. In: Journal of Geophysical Research: Oceans 96.C8,
pp. 14833–14848. ISSN: 2156-2202. DOI: 10.1029/91JC01163.
Checkley, David M. et al. (2009). “Patterns and Processes in the California Current System”. In:
Progress in Oceanography. Eastern Boundary Upwelling Ecosystems: Integrative and Comparative Approaches 83.1, pp. 49–64. ISSN: 0079-6611. DOI: 10.1016/j.pocean.2009.07.028.
Chelton, Dudley B. et al. (2011). “The Influence of Nonlinear Mesoscale Eddies on Near-Surface
Oceanic Chlorophyll”. In: Science 334.6054, pp. 328–332. DOI: 10.1126/science.1208897.
Chrysagi, Evridiki et al. (2021). “High-Resolution Simulations of Submesoscale Processes in the
Baltic Sea: The Role of Storm Events”. In: Journal of Geophysical Research: Oceans 126.3,
e2020JC016411. ISSN: 2169-9291. DOI: 10.1029/2020JC016411.
Clayton, Sophie et al. (2013). “Dispersal, Eddies, and the Diversity of Marine Phytoplankton”. In:
Limnology and Oceanography: Fluids and Environments 3.1, pp. 182–197. ISSN: 2157-3689.
DOI: 10.1215/21573689-2373515.
Clayton, Sophie et al. (2017). “Co-Existence of Distinct Ostreococcus Ecotypes at an Oceanic
Front”. In: Limnology and Oceanography 62.1, pp. 75–88. ISSN: 1939-5590. DOI: 10.1002/
lno.10373.
Cochlan, William (2008). “Nitrogen Uptake in the Southern Ocean”. In: pp. 569–596. ISBN: 978-
0-12-372522-6.
100
d’Ovidio, F. et al. (2010). “Fluid Dynamical Niches of Phytoplankton Types”. In: Proceedings of
the National Academy of Sciences 107.43, pp. 18366–18370. ISSN: 0027-8424, 1091-6490.
DOI: 10.1073/pnas.1004620107.
Damien, Pierre et al. (2023). “Modulation of Phytoplankton Uptake by Mesoscale and Submesoscale Eddies in the California Current System”. In: Geophysical Research Letters 50.16,
e2023GL104853. ISSN: 1944-8007. DOI: 10.1029/2023GL104853.
Dauhajre, Daniel P. et al. (2017). “Submesoscale Coherent Structures on the Continental Shelf”.
In: Journal of Physical Oceanography 47.12, pp. 2949–2976. ISSN: 0022-3670, 1520-0485.
DOI: 10.1175/JPO-D-16-0270.1.
Deser, Clara et al. (2010). “Sea Surface Temperature Variability: Patterns and Mechanisms”. In:
Annual Review of Marine Science 2.1, pp. 115–143. DOI: 10 . 1146 / annurev - marine -
120408-151453.
Deutsch, Curtis et al. (2021). “Biogeochemical Variability in the California Current System”. In:
Progress in Oceanography 196, p. 102565. ISSN: 0079-6611. DOI: 10 . 1016 / j . pocean .
2021.102565.
DiGiacomo, Paul M. et al. (2001). “Satellite Observations of Small Coastal Ocean Eddies in the
Southern California Bight”. In: Journal of Geophysical Research: Oceans 106.C10, pp. 22521–
22543. ISSN: 2156-2202. DOI: 10.1029/2000JC000728.
Doblin, Martina A. et al. (2016). “Drift in Ocean Currents Impacts Intergenerational Microbial Exposure to Temperature”. In: Proceedings of the National Academy of Sciences 113.20, pp. 5700–
5705. ISSN: 0027-8424, 1091-6490. DOI: 10.1073/pnas.1521093113.
Doney, Scott C. (1999). “Major Challenges Confronting Marine Biogeochemical Modeling”. In:
Global Biogeochemical Cycles 13.3, pp. 705–714. ISSN: 1944-9224. DOI: 10.1029/1999GB900039.
Dong, Changming et al. (2007). “A Numerical Study of Island Wakes in the Southern California
Bight”. In: Continental Shelf Research. Recent Developments in Physical Oceanographic Modelling: Part IV 27.9, pp. 1233–1248. ISSN: 0278-4343. DOI: 10.1016/j.csr.2007.01.016.
Dong, Changming et al. (2009). “Circulation and Multiple-Scale Variability in the Southern California Bight”. In: Progress in Oceanography 82.3, pp. 168–190. ISSN: 0079-6611. DOI: 10.
1016/j.pocean.2009.07.005.
Dugdale, R. C. et al. (1989). “New Production in the Upwelling Center at Point Conception, California: Temporal and Spatial Patterns”. In: Deep Sea Research Part A. Oceanographic Research
Papers 36.7, pp. 985–1007. ISSN: 0198-0149. DOI: 10.1016/0198-0149(89)90074-5.
Dugdale, R. C. et al. (2006). “Nutrient Controls on New Production in the Bodega Bay, California,
Coastal Upwelling Plume”. In: Deep Sea Research Part II: Topical Studies in Oceanography.
The Role of Wind-Driven Flow in Shelf Productivity 53.25, pp. 3049–3062. ISSN: 0967-0645.
DOI: 10.1016/j.dsr2.2006.07.009.
Eppley, Richard (1972). “Temperature and Phytoplankton Growth in the Sea”. In: Fishery Bulliten
70.4, pp. 1063–1085.
Ernst, Paul A. et al. (2023). “Characteristics of Submesoscale Eddy Structures within Mesoscale
Eddies in the Gulf of Mexico from 1/48◦ ECCO Estimates”. In: Frontiers in Marine Science
10. ISSN: 2296-7745.
Evans, Mary Anne et al. (2008). “Internal Wave Effects on Photosynthesis: Experiments, Theory,
and Modeling”. In: Limnology and Oceanography 53.1, pp. 339–353. ISSN: 1939-5590. DOI:
10.4319/lo.2008.53.1.0339.
101
Falkowski, Paul G. et al. (2003). “Phytoplankton and Their Role in Primary, New, and Export Production”. In: Ocean Biogeochemistry: The Role of the Ocean Carbon Cycle in Global Change.
Ed. by Michael J. R. Fasham. Berlin, Heidelberg: Springer, pp. 99–121. ISBN: 978-3-642-
55844-3. DOI: 10.1007/978-3-642-55844-3_5.
Falkowski, Paul G. et al. (2008). “The Microbial Engines That Drive Earth’s Biogeochemical Cycles”. In: Science 320.5879, pp. 1034–1039. ISSN: 0036-8075, 1095-9203. DOI: 10 . 1126 /
science.1153213.
Field, Christopher B. et al. (1998). “Primary Production of the Biosphere: Integrating Terrestrial
and Oceanic Components”. In: Science 281.5374, pp. 237–240. DOI: 10.1126/science.281.
5374.237.
Flexas, M. Mar et al. (2019). “Global Estimates of the Energy Transfer From the Wind to the
Ocean, With Emphasis on Near-Inertial Oscillations”. In: Journal of Geophysical Research:
Oceans 124.8, pp. 5723–5746. ISSN: 2169-9291. DOI: 10.1029/2018JC014453.
Follows, Michael J. et al. (2007). “Emergent Biogeography of Microbial Communities in a Model
Ocean”. In: Science 315.5820, pp. 1843–1846. ISSN: 0036-8075, 1095-9203. DOI: 10.1126/
science.1138544.
Freilich, Mara A. et al. (2019). “Decomposition of Vertical Velocity for Nutrient Transport in the
Upper Ocean”. In: Journal of Physical Oceanography 49.6, pp. 1561–1575. ISSN: 0022-3670,
1520-0485. DOI: 10.1175/JPO-D-19-0002.1.
Frischknecht, M. et al. (2018). “Origin, Transformation, and Fate: The Three-Dimensional Biological Pump in the California Current System”. In: Journal of Geophysical Research: Oceans
123.11, pp. 7939–7962. ISSN: 2169-9291. DOI: 10.1029/2018JC013934.
Fu, Weiwei et al. (2016). “Climate Change Impacts on Net Primary Production (NPP) and Export
Production(EP) Regulated by Increasing Stratification and Phytoplankton Communitystructure
in the CMIP5 Models”. In: Biogeosciences 13.18, pp. 5151–5170. ISSN: 1726-4189. DOI: 10.
5194/bg-13-5151-2016.
Garwood, Jessica C. et al. (2020). “Life in Internal Waves”. In: Oceanography 33.3, pp. 38–49.
ISSN: 10428275, 2377617X. JSTOR: 26962480.
Gaube, P. et al. (2013). “Satellite Observations of Chlorophyll, Phytoplankton Biomass, and Ekman Pumping in Nonlinear Mesoscale Eddies”. In: Journal of Geophysical Research: Oceans
118.12, pp. 6349–6370. ISSN: 2169-9291. DOI: 10.1002/2013JC009027.
Geider, Richard J. et al. (1996). “A Dynamic Model of Photoadaptation in Phytoplankton”. In:
Limnology and Oceanography 41.1, pp. 1–15. ISSN: 1939-5590. DOI: 10.4319/lo.1996.41.
1.0001.
Graham, Robert M. et al. (2013). “The Dynamical Subtropical Front”. In: Journal of Geophysical
Research: Oceans 118.10, pp. 5676–5685. ISSN: 2169-9291. DOI: 10.1002/jgrc.20408.
Guo, Lin et al. (2014). “Seasonal Dynamics of Physical and Biological Processes in the Central
California Current System: A Modeling Study”. In: Ocean Dynamics 64.8, pp. 1137–1152.
ISSN: 16167341. DOI: 10.1007/s10236-014-0721-x.
Guo, Lin et al. (2020). “Modeled Dynamics of Physical and Biological Processes in the Central
California Current System From 1993 to 2016”. In: Journal of Geophysical Research: Oceans
125.5, e2019JC015664. ISSN: 2169-9291. DOI: 10.1029/2019JC015664.
Hellweger, Ferdi L. et al. (2007). “Accounting for Intrapopulation Variability in Biogeochemical Models Using Agent-Based Methods”. In: Environmental Science & Technology 41.8,
pp. 2855–2860. ISSN: 0013-936X. DOI: 10.1021/es062046j.
102
Hellweger, Ferdi L. et al. (2016). “The Role of Ocean Currents in the Temperature Selection of
Plankton: Insights from an Individual-Based Model”. In: PLOS ONE 11.12. Ed. by Hans G.
Dam, e0167010. ISSN: 1932-6203. DOI: 10.1371/journal.pone.0167010.
Hood, Raleigh R. et al. (1990). “Surface Patterns in Temperature, Flow, Phytoplankton Biomass,
and Species Composition in the Coastal Transition Zone off Northern California”. In: Journal
of Geophysical Research: Oceans 95.C10, pp. 18081–18094. ISSN: 2156-2202. DOI: 10.1029/
JC095iC10p18081.
Hutchinson, G. E. (1961). “The Paradox of the Plankton”. In: The American Naturalist 95.882,
pp. 137–145. ISSN: 0003-0147. DOI: 10.1086/282171.
Jensen, J. L. W. V. (1906). “Sur les fonctions convexes et les inegalit ´ es entre les valeurs moyennes”. ´
In: Acta Mathematica 30, pp. 175–193. ISSN: 0001-5962, 1871-2509. DOI: 10.1007/BF02418571.
Jing, Zhiyou et al. (2020). “Submesoscale Fronts and Their Dynamical Processes Associated with
Symmetric Instability in the Northwest Pacific Subtropical Ocean”. In: Journal of Physical
Oceanography 51.1, pp. 83–100. ISSN: 0022-3670, 1520-0485. DOI: 10.1175/JPO- D- 20-
0076.1.
Jonsson, Bror F. et al. (2016). “Episodicity in Phytoplankton Dynamics in a Coastal Region”. ¨
In: Geophysical Research Letters 43.11, pp. 5821–5828. ISSN: 1944-8007. DOI: 10 . 1002 /
2016GL068683.
Kamykowski, Daniel (1974). “Possible Interactions between Phytoplankton and Semidiurnal Internal Tides”. In: Journal of Marine Research 32.1.
Kessouri, Fayc¸al et al. (2020). “Submesoscale Currents Modulate the Seasonal Cycle of Nutrients
and Productivity in the California Current System”. In: Global Biogeochemical Cycles 34.10,
e2020GB006578. ISSN: 1944-9224. DOI: 10.1029/2020GB006578.
Kida, Shinichiro et al. (2017). “A Lagrangian View of Spring Phytoplankton Blooms”. In: Journal
of Geophysical Research: Oceans 122.11, pp. 9160–9175. ISSN: 2169-9291. DOI: 10.1002/
2017JC013383.
Kling, Joshua D. et al. (2019). “Transient Exposure to Novel High Temperatures Reshapes Coastal
Phytoplankton Communities”. In: The ISME Journal. ISSN: 1751-7362, 1751-7370. DOI: 10.
1038/s41396-019-0525-6.
Kremer, Colin T. et al. (2018). “Gradual Plasticity Alters Population Dynamics in Variable Environments: Thermal Acclimation in the Green Alga Chlamydomonas Reinhartdii”. In: Proceedings of the Royal Society B: Biological Sciences 285.1870, p. 20171942. ISSN: 0962-8452,
1471-2954. DOI: 10.1098/rspb.2017.1942.
Kudela, Raphael M. et al. (2008). “New Insights into the Controls and Mechanisms of Plankton Productivity in Coastal Upwelling Waters of the Northern California Current System”. In:
Oceanography 21.4, pp. 46–59. ISSN: 1042-8275.
Lacour, Leo et al. (2015). “Phytoplankton Biomass Cycles in the North Atlantic Subpolar Gyre: A ´
Similar Mechanism for Two Different Blooms in the Labrador Sea”. In: Geophysical Research
Letters 42.13, pp. 5403–5410. ISSN: 1944-8007. DOI: 10.1002/2015GL064540.
Landry, Michael R. et al. (2009). “Lagrangian Studies of Phytoplankton Growth and Grazing Relationships in a Coastal Upwelling Ecosystem off Southern California”. In: Progress in Oceanography. Eastern Boundary Upwelling Ecosystems: Integrative and Comparative Approaches
83.1, pp. 208–216. ISSN: 0079-6611. DOI: 10.1016/j.pocean.2009.07.026.
103
Large, W. G. et al. (1994). “Oceanic Vertical Mixing: A Review and a Model with a Nonlocal
Boundary Layer Parameterization”. In: Reviews of Geophysics 32.4, pp. 363–403. ISSN: 1944-
9208. DOI: 10.1029/94RG01872.
Laufkotter, C. et al. (2015). “Drivers and Uncertainties of Future Global Marine Primary Produc- ¨
tion in Marine Ecosystem Models”. In: Biogeosciences 12.23, pp. 6955–6984. ISSN: 1726-
4189. DOI: 10.5194/bg-12-6955-2015.
Lazar, Ayah et al. (2018). “Submesoscale Turbulence over a Topographic Slope”. In: Fluids 3.2,
p. 32. ISSN: 2311-5521. DOI: 10.3390/fluids3020032.
Lehahn, Yoav et al. (2018). “A Satellite-Based Lagrangian View on Phytoplankton Dynamics”.
In: Annual Review of Marine Science 10.1, pp. 99–119. DOI: 10.1146/annurev- marine121916-063204.
Levy, M. et al. (2010). “Modifications of Gyre Circulation by Sub-Mesoscale Physics”. In: ´ Ocean
Modelling 34.1, pp. 1–15. ISSN: 1463-5003. DOI: 10.1016/j.ocemod.2010.04.001.
Levy, M. et al. (2012a). “Large-Scale Impacts of Submesoscale Dynamics on Phytoplankton: Local ´
and Remote Effects”. In: Ocean Modelling 43–44, pp. 77–93. ISSN: 1463-5003. DOI: 10.1016/
j.ocemod.2011.12.003.
Levy, Marina et al. (2012b). “Bringing Physics to Life at the Submesoscale”. In: ´ Geophysical
Research Letters 39.14. ISSN: 1944-8007. DOI: 10.1029/2012GL052756.
Levy, Marina et al. (2014). “Phytoplankton Diversity and Community Structure Affected by Oceanic ´
Dispersal and Mesoscale Turbulence”. In: Limnology and Oceanography: Fluids and Environments 4.1, pp. 67–84. ISSN: 2157-3689. DOI: 10.1215/21573689-2768549.
Levy, Marina et al. (2018). “The Role of Submesoscale Currents in Structuring Marine Ecosys- ´
tems”. In: Nature Communications 9.1, p. 4758. ISSN: 2041-1723. DOI: 10.1038/s41467-
018-07059-3.
Li, Qian P. et al. (2010). “Modeling Phytoplankton Growth Rates and Chlorophyll to Carbon Ratios
in California Coastal and Pelagic Ecosystems”. In: Journal of Geophysical Research: Biogeosciences 115.G4. ISSN: 2156-2202. DOI: 10.1029/2009JG001111.
Li, Qian P. et al. (2012). “Enhanced Nitrate Fluxes and Biological Processes at a Frontal Zone in
the Southern California Current System”. In: Journal of Plankton Research 34.9, pp. 790–801.
ISSN: 0142-7873. DOI: 10.1093/plankt/fbs006.
Lin, Hongyang et al. (2020). “Characterizing Meso- to Submesoscale Features in the South China
Sea”. In: Progress in Oceanography 188, p. 102420. ISSN: 0079-6611. DOI: 10 . 1016 / j .
pocean.2020.102420.
Liu, Guangpeng et al. (2018). “The Influence of Mesoscale and Submesoscale Circulation on Sinking Particles in the Northern Gulf of Mexico”. In: Elementa: Science of the Anthropocene 6.
Ed. by Jody W. Deming et al., p. 36. ISSN: 2325-1026. DOI: 10.1525/elementa.292.
Liu, Guangpeng et al. (2021). “Submesoscale Mixing Across the Mixed Layer in the Gulf of Mexico”. In: Frontiers in Marine Science 8. ISSN: 2296-7745. DOI: 10 . 3389 / fmars . 2021 .
615066.
Liu, Xiao et al. (2016). “Enhancement of Phytoplankton Chlorophyll by Submesoscale Frontal
Dynamics in the North Pacific Subtropical Gyre”. In: Geophysical Research Letters 43.4,
pp. 1651–1659. ISSN: 1944-8007. DOI: 10.1002/2015GL066996.
Lomas, M. W. et al. (2010). “Increased Ocean Carbon Export in the Sargasso Sea Linked to Climate
Variability Is Countered by Its Enhanced Mesopelagic Attenuation”. In: Biogeosciences 7.1,
pp. 57–70. ISSN: 1726-4189. DOI: 10.5194/bg-7-57-2010.
104
Lucas, Andrew J. et al. (2011a). “Horizontal Internal-Tide Fluxes Support Elevated Phytoplankton
Productivity over the Inner Continental Shelf”. In: Limnology and Oceanography: Fluids and
Environments 1.1, pp. 56–74. ISSN: 2157-3689. DOI: 10.1215/21573698-1258185.
Lucas, Andrew J. et al. (2011b). “The Green Ribbon: Multiscale Physical Control of Phytoplankton
Productivity and Community Structure over a Narrow Continental Shelf”. In: Limnology and
Oceanography 56.2, pp. 611–626. ISSN: 1939-5590. DOI: 10.4319/lo.2011.56.2.0611.
Lynn, Ronald J. et al. (1987). “The California Current System: The Seasonal Variability of Its
Physical Characteristics”. In: Journal of Geophysical Research: Oceans 92.C12, pp. 12947–
12966. ISSN: 2156-2202. DOI: 10.1029/JC092iC12p12947.
Lynn, Ronald J. et al. (2003). “Seasonal Renewal of the California Current: The Spring Transition
off California”. In: Journal of Geophysical Research: Oceans 108.C8. ISSN: 2156-2202. DOI:
10.1029/2003JC001787.
Mahadevan, Amala (2016). “The Impact of Submesoscale Physics on Primary Productivity of
Plankton”. In: Annual Review of Marine Science 8.1, pp. 161–184. DOI: 10.1146/annurevmarine-010814-015912.
Mahadevan, Amala et al. (2006). “An Analysis of Mechanisms for Submesoscale Vertical Motion
at Ocean Fronts”. In: Ocean Modelling 14.3, pp. 241–256. ISSN: 1463-5003. DOI: 10.1016/
j.ocemod.2006.05.006.
Mahadevan, Amala et al. (2012). “Eddy-Driven Stratification Initiates North Atlantic Spring Phytoplankton Blooms”. In: Science 337.6090, pp. 54–58. ISSN: 0036-8075, 1095-9203. DOI: 10.
1126/science.1218740.
Maheshwari, Megha et al. (2013). “An Investigation of the Southern Ocean Surface Temperature Variability Using Long-Term Optimum Interpolation SST Data”. In: ISRN Oceanography
2013, pp. 1–9. ISSN: 2090-8989. DOI: 10.5402/2013/392632.
Mangolte, Ines et al. (2022). “Plankton Community Response to Fronts: Winners and Losers”. In: `
Journal of Plankton Research 44.2, pp. 241–258. ISSN: 0142-7873. DOI: 10.1093/plankt/
fbac010.
Matthews, D. K. et al. (2009). “Velocity Observations of the California Current Derived from
Satellite Imagery”. In: Journal of Geophysical Research: Oceans 114.C8. ISSN: 2156-2202.
DOI: 10.1029/2008JC005029.
McKee, Darren C. et al. (2023). “Biophysical Dynamics at Ocean Fronts Revealed by Bio-Argo
Floats”. In: Journal of Geophysical Research: Oceans 128.3, e2022JC019226. ISSN: 2169-
9291. DOI: 10.1029/2022JC019226.
McWilliams, James C. (2016). “Submesoscale Currents in the Ocean”. In: Proceedings of the
Royal Society A: Mathematical, Physical and Engineering Sciences 472.2189, p. 20160117.
DOI: 10.1098/rspa.2016.0117.
Mitchell, Greg et al. (2017). “Photosynthetically Active Radiation (PAR) at Depths from 0 to
100 Meters, Expressed as Percentage of Surface PAR, Measured Aboard CCE LTER Process
Cruises in the California Current, 2006, 2007 and 2008.” In.
Moisan, John R. et al. (2002). “Modelling the Effect of Temperature on the Maximum Growth
Rates of Phytoplankton Populations”. In: Ecological Modelling 153.3, pp. 197–215. ISSN:
0304-3800. DOI: 10.1016/S0304-3800(02)00008-X.
Molemaker, M. Jeroen et al. (2005). “Baroclinic Instability and Loss of Balance”. In: Journal of
Physical Oceanography 35.9, pp. 1505–1517. ISSN: 0022-3670, 1520-0485. DOI: 10.1175/
JPO2770.1.
105
Morrow, Rosemary et al. (2019). “Global Observations of Fine-Scale Ocean Surface Topography
With the Surface Water and Ocean Topography (SWOT) Mission”. In: Frontiers in Marine
Science 6. ISSN: 2296-7745. DOI: 10.3389/fmars.2019.00232.
Morvan, Mathieu et al. (2019). “The Life Cycle of Submesoscale Eddies Generated by Topographic Interactions”. In: Ocean Science 15.6, pp. 1531–1543. ISSN: 1812-0784. DOI: 10 .
5194/os-15-1531-2019.
Omand, Melissa M. et al. (2017). “Using Bio-Optics to Reveal Phytoplankton Physiology from a
Wirewalker Autonomous Platform”. In: Oceanography 30.2, pp. 128–131. ISSN: 1042-8275.
JSTOR: 26201860.
Paparella, Francesco et al. (2020). “Stirring, Mixing, Growing: Microscale Processes Change
Larger Scale Phytoplankton Dynamics”. In: Frontiers in Marine Science 7. ISSN: 2296-7745.
DOI: 10.3389/fmars.2020.00654.
Patmore, Ryan D. et al. (2024). “Evaluating Existing Ocean Glider Sampling Strategies for Submesoscale Dynamics”. In: DOI: 10.1175/JTECH-D-23-0055.1.
Payandeh, A. R. et al. (2023). “The Occurrence, Variability, and Potential Drivers of Submesoscale
Eddies in the Southern California Bight Based on a Decade of High-Frequency Radar Observations”. In: Journal of Geophysical Research: Oceans 128.10, e2023JC019914. ISSN: 2169-
9291. DOI: 10.1029/2023JC019914.
Qu, Pingping et al. (2019). “Distinct Responses of the Nitrogen-Fixing Marine Cyanobacterium
Trichodesmium to a Thermally Variable Environment as a Function of Phosphorus Availability”. In: Frontiers in Microbiology 10, p. 1282. ISSN: 1664-302X. DOI: 10 . 3389 / fmicb .
2019.01282.
Quer´ e, Corinne Le et al. (2005). “Ecosystem Dynamics Based on Plankton Functional Types for ´
Global Ocean Biogeochemistry Models”. In: Global Change Biology 11.11, pp. 2016–2040.
ISSN: 1365-2486. DOI: 10.1111/j.1365-2486.2005.1004.x.
Reynolds, Richard W. et al. (1994). “Improved Global Sea Surface Temperature Analyses Using
Optimum Interpolation”. In: Journal of Climate 7.6, pp. 929–948. ISSN: 0894-8755. DOI: 10.
1175/1520-0442(1994)007<0929:IGSSTA>2.0.CO;2.
Ricour, Florian et al. (2023). “Century-Scale Carbon Sequestration Flux throughout the Ocean by
the Biological Pump”. In: Nature Geoscience 16.12, pp. 1105–1113. ISSN: 1752-0908. DOI:
10.1038/s41561-023-01318-9.
Rohr, Tyler et al. (2020a). “Eddy-Modified Iron, Light, and Phytoplankton Cell Division Rates
in the Simulated Southern Ocean”. In: Global Biogeochemical Cycles 34.6, e2019GB006380.
ISSN: 1944-9224. DOI: 10.1029/2019GB006380.
— (2020b). “The Simulated Biological Response to Southern Ocean Eddies via Biological Rate
Modification and Physical Transport”. In: Global Biogeochemical Cycles 34.6, e2019GB006385.
ISSN: 1944-9224. DOI: 10.1029/2019GB006385.
Rudnick, Daniel L. (2001). “On the Skewness of Vorticity in the Upper Ocean”. In: Geophysical
Research Letters 28.10, pp. 2045–2048. ISSN: 1944-8007. DOI: 10.1029/2000GL012265.
Rumyantseva, Anna et al. (2019). “Phytoplankton Spring Bloom Initiation: The Impact of Atmospheric Forcing and Light in the Temperate North Atlantic Ocean”. In: Progress in Oceanography 178, p. 102202. ISSN: 0079-6611. DOI: 10.1016/j.pocean.2019.102202.
Sandstrom, H. et al. (1984). “Internal Tide and Solitons on the Scotian Shelf: A Nutrient Pump at
Work”. In: Journal of Geophysical Research: Oceans 89.C4, pp. 6415–6426. ISSN: 2156-2202.
DOI: 10.1029/JC089iC04p06415.
106
Schaum, C.-Elisa et al. (2018). “Environmental Fluctuations Accelerate Molecular Evolution of
Thermal Tolerance in a Marine Diatom”. In: Nature Communications 9.1, p. 1719. ISSN: 2041-
1723. DOI: 10.1038/s41467-018-03906-5.
Sef´ erian, Roland et al. (2020). “Tracking Improvement in Simulated Marine Biogeochemistry Be- ´
tween CMIP5 and CMIP6”. In: Current Climate Change Reports 6.3, pp. 95–119. ISSN: 2198-
6061. DOI: 10.1007/s40641-020-00160-0.
Sherman, Elliot et al. (2016). “Temperature Influence on Phytoplankton Community Growth Rates”.
In: Global Biogeochemical Cycles 30.4, pp. 550–559. ISSN: 08866236. DOI: 10.1002/2015GB005272.
Shi, Weian et al. (2024). “Asymmetry of Submesoscale Instabilities in Anticyclonic and Cyclonic
Eddies”. In: Geophysical Research Letters 51.2, e2023GL106853. ISSN: 1944-8007. DOI: 10.
1029/2023GL106853.
Soccodato, Alice et al. (2016). “Estimating Planktonic Diversity through Spatial Dominance Patterns in a Model Ocean”. In: Marine Genomics 29, pp. 9–17. ISSN: 1874-7787. DOI: 10.1016/
j.margen.2016.04.015.
Strub, P. Ted et al. (1991). “The Nature of the Cold Filaments in the California Current System”.
In: Journal of Geophysical Research: Oceans 96.C8, pp. 14743–14768. ISSN: 2156-2202. DOI:
10.1029/91JC01024.
Swart, Sebastiaan et al. (2020). “Submesoscale Fronts in the Antarctic Marginal Ice Zone and
Their Response to Wind Forcing”. In: Geophysical Research Letters 47.6, e2019GL086649.
ISSN: 1944-8007. DOI: 10.1029/2019GL086649.
Taylor, Andrew G. et al. (2012). “Sharp Gradients in Phytoplankton Community Structure across
a Frontal Zone in the California Current Ecosystem”. In: Journal of Plankton Research 34.9,
pp. 778–789. ISSN: 0142-7873. DOI: 10.1093/plankt/fbs036.
Taylor, Andrew G. et al. (2018). “Phytoplankton Biomass and Size Structure across Trophic Gradients in the Southern California Current and Adjacent Ocean Ecosystems”. In: Marine Ecology
Progress Series 592, pp. 1–17. ISSN: 0171-8630, 1616-1599. DOI: 10.3354/meps12526.
Taylor, J. R. (2016). “Turbulent Mixing, Restratification, and Phytoplankton Growth at a Submesoscale Eddy”. In: Geophysical Research Letters 43.11, pp. 5784–5792. ISSN: 1944-8007. DOI:
10.1002/2016GL069106.
Taylor, John R. et al. (2020). “The Influence of Submesoscales and Vertical Mixing on the Export
of Sinking Tracers in Large-Eddy Simulations”. In: Journal of Physical Oceanography 50.5,
pp. 1319–1339. ISSN: 0022-3670, 1520-0485. DOI: 10.1175/JPO-D-19-0267.1.
Thomas, Mridul K. et al. (2012). “A Global Pattern of Thermal Adaptation in Marine Phytoplankton”. In: Science 338.6110, pp. 1085–1088. ISSN: 0036-8075, 1095-9203. DOI: 10 . 1126 /
science.1224836.
Torres, Hector et al. (2022). “Separating Energetic Internal Gravity Waves and Small-Scale Frontal
Dynamics”. In: Geophysical Research Letters. DOI: 10.1029/2021GL096249.
Torres, Hector S. et al. (2018). “Partitioning Ocean Motions Into Balanced Motions and Internal
Gravity Waves: A Modeling Study in Anticipation of Future Space Missions”. In: Journal of
Geophysical Research: Oceans 123.11, pp. 8084–8105. ISSN: 2169-9291. DOI: 10 . 1029 /
2018JC014438.
Trotta, Francesco et al. (2017). “Multi-Nest High-Resolution Model of Submesoscale Circulation
Features in the Gulf of Taranto”. In: Ocean Dynamics 67.12, pp. 1609–1625. ISSN: 1616-7228.
DOI: 10.1007/s10236-017-1110-z.
107
Van der Stocken, Tom et al. (2019). “Global-Scale Dispersal and Connectivity in Mangroves”. In:
Proceedings of the National Academy of Sciences 116.3, pp. 915–922. DOI: 10.1073/pnas.
1812470116.
Wang, Xinwei et al. (2019). “How Will the Key Marine Calcifier <I>Emiliania Huxleyi</I>
Respond to a Warmer and More Thermally Variable Ocean?” In: Biogeosciences 16.22, pp. 4393–
4409. ISSN: 1726-4189. DOI: 10.5194/bg-16-4393-2019.
Ward, Ben A. et al. (2016). “Marine Mixotrophy Increases Trophic Transfer Efficiency, Mean Organism Size, and Vertical Carbon Flux”. In: Proceedings of the National Academy of Sciences
113.11, pp. 2958–2963. DOI: 10.1073/pnas.1517118113.
Ward, Ben A. et al. (2017). “The Size Dependence of Phytoplankton Growth Rates: A Trade-Off
between Nutrient Uptake and Metabolism”. In: The American Naturalist 189.2, pp. 170–177.
ISSN: 0003-0147, 1537-5323. DOI: 10.1086/689992.
Webb, Eric A. et al. (2009). “Phenotypic and Genotypic Characterization of Multiple Strains of
the Diazotrophic Cyanobacterium, Crocosphaera Watsonii, Isolated from the Open Ocean”. In:
Environmental Microbiology 11.2, pp. 338–348. ISSN: 1462-2920. DOI: 10.1111/j.1462-
2920.2008.01771.x.
Wenegrat, Jacob O. et al. (2020). “Enhanced Mixing across the Gyre Boundary at the Gulf Stream
Front”. In: Proceedings of the National Academy of Sciences. ISSN: 0027-8424, 1091-6490.
DOI: 10.1073/pnas.2005558117.
Wilkerson, Frances P. et al. (2006). “The Phytoplankton Bloom Response to Wind Events and Upwelled Nutrients during the CoOP WEST Study”. In: Deep Sea Research Part II: Topical Studies in Oceanography. The Role of Wind-Driven Flow in Shelf Productivity 53.25, pp. 3023–
3048. ISSN: 0967-0645. DOI: 10.1016/j.dsr2.2006.07.007.
Winant, Clinton D. et al. (1997). “Seasonal Patterns of Surface Wind Stress and Heat Flux over the
Southern California Bight”. In: Journal of Geophysical Research: Oceans 102.C3, pp. 5641–
5653. ISSN: 2156-2202. DOI: 10.1029/96JC02801.
Zaiss, Jessica et al. (2021). “Impact of Lagrangian Sea Surface Temperature Variability on Southern Ocean Phytoplankton Community Growth Rates”. In: Global Biogeochemical Cycles 35.8,
e2020GB006880. ISSN: 1944-9224. DOI: 10.1029/2020GB006880.
108
Appendices
A1 Chapter 1 Supplement
A1.1 Introduction
This supplemental material contains the results of the sensitivity analyses we performed on reaction
norm width, final SST in the idealized simulations, and the magnitude of the imposed minimum
biomass, as well as the statistical analysis to determine the significance, or lack of, the differences
between the sensitivity tests and the results in the main text. Also included are figures to supplement the findings in the main text such as results for the broad shaped reaction norms and the
decrease ∆SST results for the skewed reaction norms. Tables containing SST variability analyses
as well as data on acclimation rates can also be found here.
A1.2 Impact of final temperature in the idealized simulations
To assess the sensitivity of our choice of final SST of 15 °C for the idealized simulations, we
performed 100 idealized simulations with final SSTs of 10 °C and 20 °C with the same rates and
magnitudes of temperature change as presented in the main text. Specifically, we compared the
percent difference between the individual-based model and the Eppley growth model, relative to
Eppley ([Eppley – phenotype model]/Eppley), at the timestep when SST stabilizes, as well as the
length of time it took for the community growth rates to equilibrate to steady conditions, referred
to as the memory length in the main text.
The final SST of the idealized profiles did not impact the results of our study. The offset between the Eppley growth model and the phenotype model were not different, statistically (95% CI,
109
see section S.7 for description of statistical analyses) for all the simulations (Figure S2). Additionally, because the results are presented in terms of generations rather than absolute time, the length
of the memory effect is also statistically not different (95% CI) for all simulations (Figure S3).
A1.3 Impact of reaction norm width in idealized simulations
The width of the reaction norm, or thermal niche, had a varying, but predictable, impact on the
percent difference from the Eppley growth model, relative to the Eppley model (calculated as in
S.1.). To test the impact of the reaction norm width, we ran 100 simulations of the same idealized
SST profiles as in the main text with both narrower and wider reaction norms for the individuals.
Making the thermal niche narrower (10.5 °C) relative to the simulations in the main text (thermal
niche = 14 °C) did not have a significant (95% CI) impact on the percent difference between the
individual-based model and the Eppley growth model, for either shape of reaction norm (Figure
S4, left column). Increasing the thermal niche from 14 °C to 20.5 °C (Figure S4, right column)
also did not have a significant impact on the offset from the Eppley growth model for skewed
reaction norms (95% CI) but did for the broad reaction norms (95% CI). A wider reaction norm
for the broad reaction norms decreased the percent difference from the Eppley growth model by an
average of 10.1% with simulations in which ∆SST changed over 7 days experiencing the largest
decreases (up to 29.4%).
The width of the thermal niche in conjunction with the magnitude of SST change impacted the
time lag in the community response, i.e., the memory length. For small (2-3 °C) and large (8-9
°C) SST changes, wider thermal niches produced shorter memory effects (95% CI) by an average
of 2.1 generations ±3.4 generations (1σ) for the broad reaction norms and 0.6 generations ±1.6
generations (1σ) for the skewed reaction norms. Conversely, wider thermal niches that experienced
moderate SST changes (4-7 °C) had longer memory effects than the default thermal niches by 1.8
generations ±4.0 generations (1σ) for the broad and 0.5 generations ±1.0 generations (1σ) for the
skewed reaction norms. This was seen across both sets of simulations for the broad reaction norms
(Figure S5, bottom row). For skewed reaction norms, decreasing the thermal niche width relative
110
to the default width did not have a significant impact on the memory length (95% CI) (Figure S5,
top row). Regardless of reaction norm width, the overall relationship between memory length and
the rate and magnitude of SST change was not different from the simulations in the main text and
did not change our results or conclusions.
For broad shaped reaction norms, wider reaction norms meant that individuals were able to continue to grow over a larger span of SSTs on either side of their optimum growth temperature (Topt)
compared to individuals with narrower reaction norms. When temperature changes were small (2-
3 °C), the biomass weighted community growth rate was able to better track small changes in SST
because the SST did not go outside of the thermal niche. Similarly, for large SST changes (8-9
°C), the community was able to respond to the SST changes more quickly than a community with
a narrower reaction norm because more individuals were able to grow at the final SST. When SST
changes were more moderate (4-6 °C), the individuals in the original environment could continue
to grow over a larger range of temperatures past their Topt which meant that those best suited for the
new environment had more biomass to overcome before making a significant contribution to the
biomass-weighted community growth rates compared to individuals in a community with narrow
reaction norms. The memory length for the skewed reaction norms was less affected by the width
of the reaction norm due to the asymmetry of the reaction norm shape. By increasing the width the
reaction norm, but keeping the maximum growth rate and Topt the same, the part of the reaction
norm that was extended corresponded to relatively low growth rates. So even though individuals
could grow at a larger temperature range, that growth did not have a large impact on the memory
length.
The impact of reaction norm width on the difference from the Eppley growth model and the
length of the memory effect did not change any of the conclusions of the manuscript. Across all the
simulations, larger and faster SST changes resulted in the largest offsets between the phenotype
model and the Eppley growth model and moderate SST changes induced the longest memory
effects.
111
A1.4 Model sensitivity to minimum biomass parameter
In the main text, we imposed a minimum biomass of 0.001 mmol C m-3 such that no individual
was allowed to go extinct, akin to the “everything is everywhere” principle (Hutchinson, 1961). To
test the sensitivity of our results to this parameter, we ran 100 simulations with the same idealized
SST profiles with a minimum biomass of 0.0001 mmol C m-3, an order of magnitude smaller.
For both the skewed shaped and broad shaped reaction norms, lower minimum biomass generally
increased both the offset from Q10 simulated growth rates (95% CI, Figure S5) and the memory
length (95% CI, Figure S6). The difference from Q10 increased by an average of 1.5% ±8.6%
(1σ) for the broad and 2.6% ±8.9% (1σ) for skewed shaped reaction norms, but ranged as high as
31.7% (broad) and 27.3% (skewed). For small ∆SSTs (2-3 °C), lower minimum biomass slightly
decreased the difference between Q10 simulated growth but as ∆SSTs increased, so did the offset.
Memory lengths increased by an average of 4.0 ±4.1 (1σ) generations for the broad reaction norms
and 3.0 ±3 (1σ) generations for the skewed reaction norms, but ranged as high as 12.6 generations
(broad) and 10.6 generations (skewed) longer for the smaller minimum biomass. A lower minimum
biomass meant that individuals with the minimum biomass contributed less to the overall biomassweighted community growth rate, resulting in lower growth rates and larger departures from Q10.
This also meant that those individuals best suited for the new environment started growing with
lower biomass and thus took longer to overcome the previously accumulated biomass from the
initial conditions which resulted in longer memory lengths. As such, the results presented in the
main text are a conservative estimate of the difference from Q10 and memory length.
The overall patterns remained the same between both minimum biomass simulations. The
direction of the ∆SST change did not impact the memory length for the broad reaction norms
whereas decreasing ∆SSTs yielded longer memory lengths for the skewed reaction norms for both
sets of simulations. In both sets of simulations, the moderate ∆SSTs resulted in the longest memory
lengths.
112
A1.5 Comparison of Ecosystem Model Choice
We compared the community growth rates from several different models to ensure that the results
we found were not the result of our choice of model. We found that all models showed similar
responses in community growth rate. Below is a description of each of the models used in this
comparison.
The biomass of each individual (Pi
, in mmol C m−3
) was calculated as
dPi
dt
= µi ∗Pi −loss (A1.1)
where µi (day−1
) is the individual growth rate at time t. Here we investigated different formulations for the loss term.
A1.5.1 Linear Mortality
We started with simple linear mortality, where loss scales linearly with biomass, similar to (Moisan
et al., 2002).
dPi
dt
= µi ∗Pi −m∗Pi (A1.2)
We found that, the mortality had to be set to unrealistic values (approx. equal to Q10 values)
in order to keep biomass from exponentially increasing. However, this model does still show a dip
in community growth rates with changes in SST that is described in the main text.
A1.5.2 Quadratic Mortality (used in the main text)
A more common approach is to represent loss as a quadratic mortality:
dPi
dt
= µi ∗Pi −m∗P
2
i
(A1.3)
113
Simulating phytoplankton loss as quadratic mortality showed the same dip in community
growth rates as SSTs begin to as described in the main text. The overall magnitude of the loss
term is consistent with the other models also.
A1.5.3 Simple Ecosystem
We also tested a more complex ecosystem model with linear mortality and loss due to grazing.
dPi
dt
= µi ∗Pi −g ∗
Pi
P
∗Z ∗Pi (A1.4)
where g is the temperature dependent grazing (m3 mmol C−1 day−1
) and Z is the total zooplankton biomass (mmol C m−3
). To keep our phytoplankton and zooplankton growth internally
consistent, we simultaneously solve for the change in total phytoplankton biomass (P) and zooplankton biomass (Z) over time (where P = ∑Pi for individuals whose biomass is greater than the
minimum) using the following equations:
dP
dt
= λ ∗P−g ∗Z ∗P (A1.5)
dZ
dt
= 0.3 ∗ g ∗Z ∗P−mz ∗Z (A1.6)
where λ (day−1
) is the biomass weighted community growth rate from all Pi > minimum
biomass, 0.3 is the zooplankton efficiency, and mz
is the zooplankton mortality rate (day−1
). Resolving for total P instead of using the sum of the individual biomasses allowed us to avoid issues
with resetting low biomass individuals to the minimum biomass which constantly adds biomass
to the system. This resulted in predator-prey oscillations (Figure S7) but also showed the dip in
community growth rates as SSTs began to change.
114
A1.5.4 Constant grazing
This model followed the same equations outlined for the Simple Ecosystem model above, but
instead of solving for how zooplankton biomass changes over time, we calculate Z for each time
step as:
Z = −0.0187 ∗ λ +5 ∗
λ
a
(A1.7)
where λ (day−1
) is the community growth rate defined in the main text (Equation 5) and a
(day−1
) is the growth rate from the Q10 parameter (Equation 2). This formulation provided a
relatively constant grazing pressure which prevented predator-prey oscillations. As seen with the
other formulations, this resulted in a decrease in community growth rates as SSTs change.
A1.6 Statistics Calculations for Sensitivity Tests
To calculate the potential significance of results from the sensitivity tests, we performed Type II
linear regression and tested the significance of the slope against a value of 1. The regression was
performed using the lsqfitma function in Matlab made available from the Monterey Bay Aquarium
Research Institute (https://www.mbari.org/index-of-downloadable-files/). This provided a slope and the uncertainty on that slope. Using these data, we then calculated the Z test
statistic as:
Z =
x− µ
σ
q
1
N
(A1.8)
where x is the slope to test against, here set to one, µ is the slope from the Type II regression, σ
is the standard deviation on the slope, and N is the number of independent tests to find µ, which is
one for the lsafitma regression. Once Z is calculated, we compare this to the standard score based
on a 95% confidence interval which corresponds to a standard score of ±1.96. If Z is outside of this
range, we reject the null hypothesis that the slope, µ is equal to one. Otherwise, we fail to reject
the null hypothesis.
115
A1.7 Nitrate limitation model
Phytoplankton growth is often co-limited by multiple factors. To test the impact of temperature
and nutrient co-limitation on the results presented in this study, we conducted an additional set of
model simulations. Specifically, we test the impact of including both temperature and nitrate limitation on growth. To our knowledge, there is no source of high-resolution mixed-layer nitrate data
along a Lagrangian trajectory in the Southern Ocean. While Bio-ARGO floats measure nitrate,
these floats do not accurately quantify the environmental variability experienced by phytoplankton
as they rest at depth and so are transported by deep, rather than surface, currents. Here we leverage
the predictable relationship between nitrate and temperature (Figure S9) to generate a companion
high-resolution nitrate dataset to the drifter SST dataset. Specifically, nitrate concentrations are
estimated based on the observed relationship between SST and nitrate in the upper 50m from profiling Bio-ARGO floats in the Southern Ocean (8805 data points, R2=0.95, p¡0.001). Data was
Bio-ARGO accessed through the Southern Ocean Carbon and Climate Observations and Modeling (SOCCOM) Project website (https://www.mbari.org/science/upper-ocean-systems/
chemical-sensor-group/floatviz/ Accessed: 03/08/2021). We modify the model equations
(Eq. 1 in the main text) to include nitrate limitation as:
µ(T) = α e
bT
1−
2(T −Topt)
w
2
∗Nlim (A1.9)
where nitrate limitation (Nlim) is calculated as:
Nlim =
Cnit
0.5+Cnit
(A1.10)
and Cnit is the nitrogen concentration.
When temperatures are less than 14 °C, the inclusion of both temperature and nitrate limitation
produced similar results to the model simulations with temperature only (Figure S11, Example
Trajectory 1). This is because nitrogen limitation is moderate (Nlim > 0.82) for these temperatures
and acts to decrease growth rates only slightly. For regions warmer than 14°C, nitrate limitation
116
can become substantial. Due to the non-linear relationship between nutrient concentration and
growth limitation, the onset of severe nutrient limitation is rapid (Figure S11, Example Trajectory
2). Transition from mild nutrient limitation to severe nutrient limitation results in a rapid drop in
growth rates, as expected. However, once the populations are fairly uniformly limited (i.e. stay
in waters with low nitrate) the model dynamics are once again similar to temperature limitation
alone albeit with much lower growth rates and therefore much longer memory lengths (Figure
S11, Example Trajectory 3).
117
Table A1.1: Results of SSTmax variability analysis
118
Table A1.2: Estimates of acclimation times (lower bounds) for the polar diatom F. cylindrus for
lab-cultures initiated at 2.5 °C and then exposed to first to either cooling or warming of ≈1.5 °C
and 3 °C for the first acclimation, and then subsequent step increases or decreases in temperature,
as part of the preparation for a thermal reaction norm experiment. The estimates growth rates were
used in conjunction with the change in temperature to compute the range of acclimation times
presented. For further information contact philip.boyd@utas.edu.au.
Growth
Temperature (oC)
Condition
Temperature
Change (oC)
Growth rate
(day-1) after 4
generations
Acclimation
rate (oC day-1)
-1.4 Shift down -1.3 0.2 0.26
-0.1 Initial 0.24
-0.1 Shift down -1.5 0.24 0.36
1.4 Initial 0.26
1.4 Shift down -1.4 0.26 0.36
2.8 Initial 0.26
4.1 Shift up 1.3 0.27 0.36
4.1 Initial 0.27
5.3 Shift up 1.2 0.29 0.35
5.3 Initial 0.29
6.5 Shift up 1.2 0.22 0.28
119
Figure A1.1: Box plots of the percent of the SST variability in the drifter trajectory that is accounted for by the smoothed spline. Each of the 2,190 90-day drifter and spline trajectories was
broken up into windows in 1-day increments from 1 to 90 days. The standard deviation of the
drifter trajectory is the sum of the standard deviation of the smoothed spline plus some noise term.
From this, the variability accounted for by the spline for each window, for each trajectory was
recorded with the results shown. As expected, over longer window lengths the spline accounts for
higher percentage of the overall variability.
120
Figure A1.2: The impact of final SST on percent difference between the individual based model
and the Eppley growth model, relative to Eppley growth ([Eppley – phenotype]/Eppley). The
results from the simulations in the main text are compared to simulations with final SSTs of 10 °C
(a,c) and 20 °C (b,d) for both the skewed (top row) and broad (bottom row) shaped reaction norms.
Open data points represent decreasing ∆SSTs and filled in data are increasing ∆SSTs. The black
line indicates the 1-1 line. There is no statistical difference between simulations with differing
final SSTs (95% CI).
121
Figure A1.3: The impact of final SST on memory length. The results from the simulations in the
main text are compared to simulations with final SSTs of 10 °C (a,c) and 20 °C (b,d) for both
the skewed (top row) and broad (bottom row) shaped reaction norms. Open data points represent
decreasing ∆SSTs and filled in data are increasing ∆SSTs. The black line is the 1-1 line. There is
no statistical difference between simulations with differing final SSTs (95% CI).
122
Figure A1.4: The impact of reaction norm width on percent difference from the Eppley growth
model. The results from the simulations in the main text are compared to simulations with narrower
(a,c) and wider (b,d) reaction norms for both the skewed (top row) and broad (bottom row) shaped
reaction norms. The black line is the 1-1 line. Closed data points represent increasing ∆SSTs and
open circles represent decreasing ∆SSTs. For broad reaction norms, increasing the reaction norm
width increases the difference between the phenotype model and the Eppley growth model. There
was no significant difference between the simulations for skewed reaction norms.
123
Figure A1.5: The impact of reaction norm width on memory length. The results from the simulations in the main text are compared to simulations with narrower (left column) and wider (right
column) reaction norm widths for both the skewed (top row) and broad (bottom row) shaped reaction norms. The black line is the 1-1 line. Closed data points represent increasing ∆SSTs and
open circles represent decreasing ∆SSTs. Broad reaction norms are most affected by reaction norm
width, but increasing the reaction norm width had a significant (95% CI) impact on the memory
length for the skewed reaction norms as well. For small ∆SSTs (2-3o°C), narrower reaction norms
have longer memory lengths. When ∆SSTs are large (8-9 °C), the memory length is shorter for
communities with more narrow reaction norms. For moderate ∆SST changes, (4-6 °C), the width
of the reaction norm has minimal impact on the memory length. The memory lengths associated
with moderate ∆SSTs are typically the longest memory lengths, just as in the main text.
124
Figure A1.6: The impact of minimum biomass on deviation from Q10. The results from the
simulations in the main text (x-axis) are compared to simulations with an order of magnitude
smaller minimum biomass for both skewed (left) and broad (right) shaped reaction norms. The
black line is the 1-1 line. Filled in data points represent increasing ∆SSTs and open data points
are decreasing ∆SSTs. The minimum biomass impact is significant at the 95% CI with an average
increase in offset from Q10, for both reaction norm shapes.
Figure A1.7: The impact of minimum biomass on memory length. The results from the simulations
in the main text (x-axis) are compared to simulations with an order of magnitude smaller minimum
biomass for ‘both skewed (left) and broad (right) shaped reaction norms. The black line is the 1-1
line. Filled in data points represent increasing ∆SSTs and open data points are decreasing ∆SSTs.
The minimum biomass impact is significant at the 95% CI with an average increase in memory
length for both reaction norm shapes. However, the pattern of moderate ∆SSTs exhibiting the
longest memory effects were robust across all simulations.
125
Figure A1.8: Comparison between different ecosystem model results for community growth for an
idealized simulation with an increase of 4 °C over 21 days. For community growth rates, all models
show similar qualitative results indicating a decrease in growth rate over the transient conditions
culminating in a growth rate minimum when SSTs stabilize.
126
Figure A1.9: Nitrate-temperature relationship. Relationship between temperature and nitrate in
the Southern Ocean as measured by Bio-ARGO floats. Data was collected from the Southern
Ocean Carbon and Climate Observations and Modeling (SOCCOM) Project funded by National
Science Foundation, Division of Polar Programs (NSF PLR -1425989), supplemented by NOAA
and NASA.
127
Figure A1.10: Nitrate limitation versus temperature. Due to the non-linear relationship, for temperatures greater than 14 °C, nitrate concentrations minimally impact growth rates. Once SST goes
above 14 °C, the limitation rapidly decreases to the minimum which results in large growth rate
decreases.
128
Figure A1.11: Examples of the impact of nitrate limitation on community growth rates. When
SSTs are below 14 °C, the growth rate with nitrate limitation and that without limitation are nearly
identical (Example Trajectory 1). As growth rates approach 14 °C, the growth rate dynamics become more variable but follow the same pattern as that with temperature only limitation (Example
Trajectory 2). When SSTs exceed 14 °C, the growth rate dynamics are muted due to lower overall
growth rates but still mimic those of the temperature limited only growth rates (Example Trajectory 3).
129
Figure A1.12: Probability density functions of the absolute value of the maximum change in SST
over 7, 21, 45, and 90 days for the drifter trajectories (blue) and the smoothed splines of the
trajectory SSTs (red).
130
Figure A1.13: The standard deviation (1σ) as a function of ∆SSTmax over different ∆tmax windows.
∆SSTmax drives the variability across the ∆tmax window lengths.
131
Figure A1.14: Results from SST variability analyses for the drifter (left) and satellite trajectories
(right) showing most common SST changes for each time window. Data are presented as total
percent of data that fall within that ∆SST bin for the window length. Each row sums to 100%.
Although the magnitudes of variability are similar, the nature of that variability is different with the
Lagrangian reference frame experiencing more variability consistent with longer memory effects.
132
Figure A1.15: Daily rates of SST change for drifter trajectories. The rates of change were calculated as the range of the recorded SST values over a 1-day moving window for a total of n =
781,749 data points for 197,100 days.
133
Figure A1.16: Community growth rates for each of the 100 simulations (grey lines) for an increase
of 4 °C over 7, 21, 45, and 90 days for skewed shaped reaction norms. The black line is the
Q10 simulated community growth rate and the blue line is the SST profile for the simulation.
The locations marked a1 and a2 in the 21-day panel represent the timesteps used to calculate
the percent change in growth rate associated with transient SSTs as shown in Figures 2b. This
metric was calculated as (a1 – a2)*100/a1. The locations marked b1 and b2 in the 45-day panel
represent the timesteps used to calculate the percent difference in growth rates between the Q10
parameterized growth and the phenotype model as shown in Figure 2c, S11. This metric was
calculated as (b1 – b2)*100/b1. The locations marked c1 and c2 in the 90-day panel point to the
timesteps used to calculate the memory length. The dashed grey lines represent ±5% of the final,
stabilized community growth rate which was used as the threshold for the memory effect which
was defined as the time in days between c1 when SSTs stabilize and c2 when the community
growth rate crosses the threshold.
134
Figure A1.17: Full results for the change in growth rate in the idealized simulations for the skewed
shaped reaction norms (top row) and the broad shaped reaction norms (bottom row) under both
increasing SSTs (left column) and decreasing SSTs (right column).
135
Figure A1.18: Full results for the percent difference from the Eppley growth approximation at the
time step when SSTs stabilize in the idealized simulations for the skewed shaped reaction norms
(top row) and the broad shaped reaction norms (bottom row) under both increasing SSTs (left
column) and decreasing SSTs (right column).
136
Figure A1.19: Full results for the length of the memory effect in the idealized simulations for the
skewed shaped reaction norms (top row) and the broad shaped reaction norms (bottom row) under
both increasing SSTs (left column) and decreasing SSTs (right column).
137
Figure A1.20: The 90-day average percent difference between community growth rates determined
via the Q10 method and the phenotype-based model versus standard deviation (1σ) of SST over
the 90-day trajectory. Drifter data are represented by circles colored according to their mean SST.
Black diamonds represent the first 90 days of the idealized trajectories; filled diamonds are the
idealized trajectories for which SSTs increase and open black diamonds are idealized trajectories
with decreasing SSTs. Pink triangles represent the two example trajectories from Figure 1 in the
main text.
138
Figure A1.21: Comparison of mean community growth rate over the entire 90-day trajectory for
the real trajectories and their spline simulations for skewed (left) and broad (right) shaped reaction
norms. With each reaction norm shape, smoothing the small-scale noise did not impact overall
biomass-weighted community growth rates (95% CI, t-test) further supporting that small-scale
noise does not induce a memory effect.
139
Figure A1.22: Example of impact of acclimation. (Top) Example idealized SST trajectory of
changing 2 °C in 7 days with acclimation rates of 0.2 °C day-1 and 0.3 °C day-1. Other acclimation
rates not shown as they plot along the No Acclimation line because those rates are faster than the
rate of change. (Middle) Community growth rates for skewed reaction norms over each of the
simulations for the no acclimation simulations (grey lines) and the simulations with an acclimation
rate of 0.2 °C day-1 (green lines). Dashed lines represent the thresholds used to calculate the
memory length. (Bottom) Same as the middle panel but for broad reaction norms.
140
Figure A1.23: Impact of acclimation on memory length on the skewed reaction norms (top row)
and the broad reaction norms (bottom row) in both the increasing (left column) and decreasing
(right column) ∆SST directions. Dashed lines represent the memory lengths calculated for the
simulations that did not incorporate acclimation. When acclimation rates are greater than or equal
to the rate of SST change, there is no difference between the simulations that incorporated acclimation and those that did not.
141
Figure A1.24: Comparison of the final SST for the drifter and the satellite data. The data from both
sources represent the same location in space and time so the data should be similar and in fact, are
not statistically different from one another (ttest, 95% CI). The grey line represents the 1-1 line.
142
Figure A1.25: Same as Figure 4 in the main text but for broad shaped reaction norms. The impact
of Lagrangian and Eulerian variability on community composition. Here we plot the difference
between the Topt of the most abundant phenotype at the end of each 90-day trajectory and the
final SST for the drifter trajectory (x-axis) and the satellite data (y-axis). The final SSTs for the
drifter and satellite data are statistically identical (t-test, 95% CI). Therefore, deviations from the
1:1 line demonstrate the impact of a Lagrangian versus Eulerian reference frame on community
composition.
143
Figure A1.26: The impact of SST variability on community composition. (Top) An example 90-
day drifter trajectory and the satellite SST data for the final location of the drifter over the same
90 days shown as solid lines. The dashed lines are the Topt of the most abundant phenotype at
each timestep. (Bottom) The biomass of each phenotype with a skewed shaped reaction norms at
day 90 for the satellite and drifter trajectories. In this example, the offset between the final SST is
-0.60 °C for the drifter and -0.09 °C for the satellite data. The difference in the magnitude of the
offset between the two data sets represents the difference in the variability of the SSTs. However,
in this example, the satellite SSTs stay relatively constant whereas the drifter SSTs experience a
rapid increase of 3.5 °C in 4 days beginning Sept. 29. Because the drifter SSTs remain relatively
constant through the end of the 90 days, the community is able to adjust to the new environment
before the end of the simulation which results a community Topt that reflects the SSTs at day 90
for both the satellite and the drifters.
144
Figure A1.27: The impact of SST variability on community composition. (Top) An example 90-
day drifter trajectory and the satellite SST data for the final location of the drifter over the same 90
days shown as solid lines. The dashed lines are the Topt of the most abundant phenotype at each
timestep. (Bottom) The biomass of each phenotype with a skewed shaped reaction norms at day
90 for the satellite and drifter trajectories. In this example, the offset between the final SST is -0.23
°C for the drifter and -3.1 °C for the satellite data. The difference in the magnitude of the offset
between the two data sets represents the difference in the variability of the SSTs. Here, the drifter
SSTs gradually increase over the 90 days which allows the community to continuously track the
changes in SST whereas the satellite SSTs are relatively stable and then rapidly decrease from 17.7
°C on March 10 to 13.8 °C on March 17. Due to the long memory effect associated with this rate
and magnitude of change, the community was not able to track the SST change which resulted in
a large offset between the final SST and the Topt of the most abundant phenotype at day 90.
145
Figure A1.28: Percent reduction in community growth rate from the phenotype model from the
Eppley growth model for the drifters, the satellite derived trajectories, and the satellite point data.
In the Lagrangian reference frame (both sets of trajectories) the community growth rate from the
phenotype model is lower than in the Eulerian reference frame (satellite point data) as a result of
the different SST variability encountered in each reference frame.
146
Figure A1.29: Full drifter results from Figure 9 in the main text. This figure includes drifter
trajectories that do not overlap in time with the satellite data. The spatial patterns are the same as
in the main text.
147
A2 Chapter 2 Supplement
A2.1 Results sensitivity to initial particle depth
Releasing the Lagrangian particles at different depths does not impact the interpretation of the results. The most significant difference between the simulations initialized at 100m and 50m are that
the particles released at 50m spend more time in the mixed layer. However, particles initialized at
50m encountered submesoscale features at the same relative proportions as the the particles initialized at 100m. Additionally, due to increased time in the mixed layer, the particles initialized at
50m encounter submesoscale features at a lower rate compared to the particles that were initialized
at 100m, despite the overall number of encounters being higher for the particles released at 50m.
0 0.01 0.02 0.03 0.04 0.05
Occurrence Rate by Total Area
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0.04
0.045
0.05
Encounter Rate by drifters in mixed layer
cyclonic
anticyclonic
front
50m Release
0 0.01 0.02 0.03 0.04 0.05
Occurrence Rate by Total Area
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0.04
0.045
0.05
Encounter Rate by drifters in mixed layer
cyclonic
anticyclonic
front
100m Release
0 0.01 0.02 0.03 0.04 0.05
Occurrence Rate by Total Area
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0.04
0.045
0.05
Encounter Rate by drifters in mixed layer
cyclonic
anticyclonic
front
10m Release
Figure A2.1: Spatial Occurrence rate of submesoscale features versus particle encounter rate for Lagrangian particles initialized at 100m, 50m, and 10m. Closed data points represent data from the onshore
regions. Open data points are from the offshore regions. Triangles are for the Southern regions and squares
represent the northern regions.
A2.2 Vertical Velocity Filter
Vertical velocities in the CCS are strongly influenced by large-scale dynamics such as internal
gravity waves Torres et al., 2022. These features have similar vertical velocities and timescales as
submesoscale features. To isolate the vertical velocities associated with submesoscale features in
our study, a high-pass, linear filter was applied following Torres et al., 2022.
148
Because submesoscale features are largely confined to the mixed layer, only velocities within
the mixed layer have the filter applied. For each time step, at each grid point, the depth of the
mixed layer was defined as the deepest mixing point over a 24-hour period according to the KProfile Parameterization (Large et al., 1994).
The equation that was applied over the mixed layer was:
W′ = W −m∗Z +a. (A2.1)
W was the raw unfiltered vertical velocities. Z represents the grid point depths. The variables m
and a were the slope and the intercept of the high-pass filter, respectively. The slope was calculated
as
m =
WMLD −Wo
ZMLD
. (A2.2)
WMLD was the vertical velocity at the base of the mixed layer, Wo represents the surface vertical
velocity, and ZMLD is the depth of the mixed layer. The intercept is simply the vertical velocity as
the surface.
a = Wo (A2.3)
149
Figure A2.2: The impact of filtering out large-scale vertical velocities to isolate submesoscale features.
(left) Raw, unfiltered vertical velocity at a depth of 9.67m. (middle) Filtered vertical velocities at 9.67m
depth. Submesoscale filaments and eddies are revealed. (right) Rossby number at 9.67m. Submesoscale
features are identified by absolute values >0.75. The locations of submesoscale features aligns with the
vertical velocities revealed after filtering the raw data.
A2.3 Figures
150
Figure A2.3: Standard deviation (2σ) of diffusivity over the whole simulation for a depth of 39.725 m.
151
1 2 3 4 5 6 7 8 9 10
Average size per feature [km2
]
0
50
100
150
200
250
300
Count
South Onshore
1 2 3 4 5 6 7 8 9 10
Average size per feature [km2]
0
50
100
150
200
250
300
Count
North Onshore
1 2 3 4 5 6 7 8 9 10
Average size per feature [km2]
0
20
40
60
80
100
120
140
160
Count
North Offshore
cyclonic
anticyclonic
front
1 2 3 4 5 6 7 8 9 10
Average size per feature [km2]
0
20
40
60
80
100
120
140
160
180
Count
South Offshore
Figure A2.4: The size of submesoscale features varied significantly among regions and feature types.
Cyclonic eddies were consistently the largest, averaging 2.36 ±0.73 km2
(2σ), followed by fronts (1.32
±0.43 km2
), and anticyclonic eddies (1.15 ±0.42 km2
) (Figure S2). The smaller size for anticyclonic eddies
aligns with regional observations (Payandeh et al., 2023) and is likely attributed to their susceptibility to
instabilities that limits their growth (Shi et al., 2024). Notably, features in the southern regions were, on
average, 29% larger than their northern region counter parts. Despite some overlap in size distributions,
a two-sided t-test confirmed significant differences in mean feature size among regions and feature types
(p<0.05).
152
-80 -60 -40 -20 0 20 40 60
Vertical Velocity [m/day]
0
0.02
0.04
0.06
0.08
0.10
0.12
0.14
0.16
0.18
Probability
Eulerian Mean Vertical Velocities
Cyclonic
Anticyclonic
Front
Background
-80 -60 -40 -20 0 20 40 60
Vertical Velocity [m/day]
0
0.02
0.04
0.06
0.08
0.10
0.12
Probability
Lagrangian Mean Vertical Velocities
Cyclonic
Anticyclonic
Front
Background
Figure A2.5: Distribution of submesoscale feature and background vertical velocities. (Left) Eulerian
distribution (Right) Lagrangian distribution. The mean vertical velocities were statistically different form
each other (two-sided t-test, p<0.05) and from the background average.
0 0.2 0.4 0.6 0.8 1
Horizontal Velocity [m/s]
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.10
Probability
Eulerian Mean Horizontal Velocities
Background
Cyclonic
Front
Anticyclonic
0 0.2 0.4 0.6 0.8 1
Horizontal Velocity [m/s]
0
0.005
0.010
0.015
0.020
0.025
0.030
0.035
0.040
0.045
0.050
Probability
Lagrangian Mean Horizontal Velocities
Background
Cyclonic
Front
Anticyclonic
Figure A2.6: Distribution of submesoscale feature and background horizontal velocities(Left) Eulerian
distribution (Right) Lagrangian distribution
153
-80 -60 -40 -20 0 20 40 60
Vertical Velocity [m/day]
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
Probability
North Offshore Eulerian Mean
Cyclonic
Anti
Front
Background
-80 -60 -40 -20 0 20 40 60
Vertical Velocity [m/day]
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
Probability
North Onshore Eulerian Mean
Cyclonic
Anti
Front
Background
-80 -60 -40 -20 0 20 40 60
Vertical Velocity [m/day]
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
Probability
South Offshore Eulerian Mean
Cyclonic
Anti
Front
Background
-80 -60 -40 -20 0 20 40 60
Vertical Velocity [m/day]
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
Probability
South Onshore Eulerian Mean
Cyclonic
Anti
Front
Background
Figure A2.7: Distribution of mean Eulerian vertical velocities for submesoscale features and the background. Data represent the filtered vertical velocities to remove large-scale contributions.
154
-80 -60 -40 -20 0 20 40 60
Vertical Velocity [m/day]
0
0.05
0.1
0.15
0.2
0.25
0.3
Probability
North Offshore Lagrangian Velocities
Background
Cyclonic
Anti
Front
-80 -60 -40 -20 0 20 40 60
Vertical Velocity [m/day]
0
0.05
0.1
0.15
0.2
0.25
0.3
Probability
North Onshore Lagrangian Velocities
Background
Cyclonic
Anti
Front
-80 -60 -40 -20 0 20 40 60
Vertical Velocity [m/day]
0
0.05
0.1
0.15
0.2
0.25
0.3
Probability
South Offshore Lagrangian Velocities
Background
Cyclonic
Anti
Front
-80 -60 -40 -20 0 20 40 60
Vertical Velocity [m/day]
0
0.05
0.1
0.15
0.2
0.25
0.3
Probability
South Onshore Lagrangian Velocities
Background
Cyclonic
Anti
Front
Figure A2.8: Distribution of Lagrangian vertical velocities for submesoscale features and the background.
Data represent the filtered vertical velocities to remove large-scale contributions.
155
A3 Chapter 3 Supplement
A3.1 Figures and Tables
Figure A3.1: Schematic of the quota model and the 1D, Lagrangian water column model
156
Table A3.1: Biology Model Parameters
Type Parameter Symbol Value Units Reference
Temperature limitation gammaT unitless Thomas et al., 2012
tau 0.8 unitless Thomas et al., 2012
A -4000 Kelvin Thomas et al., 2012
To 298.15 Kelvin Thomas et al., 2012
Cell size Vcell 4.20E+06 um3
Cellular Quota
carbon Qc 6.04E-05 mmol C cell-1
Nitrogen Uptake
maximum uptake
VNO3max 5.02E-07 mmol N cell-1 d-1 Ward et al., 2012
VNH4max 2.56E-07 mmol N cell-1 d-1 Ward et al., 2012
half saturation constant
kNO3 0.45 mmol NO3 m-3 Follows et al., 2007
kNH4 0.25 mmol NH4 m-3
Ammonia inhibition psi 4.6 (mmol N m-3)
-1 Ward and Follows, 2016
Photosynthesis
maximum growth rate mumax 6 day-1
basal respiration basres 0.001 day-1
initial photosynthesisirradiance slope alpha 0.05 mmol C m2
(mgChl)-1 W-1 day-1 modified from Li et al., 2010
photoinhibition beta 5.00E-05 m
2
W-1 day-1 modified from Li et al., 2010
Chl Synthesis
Chl: N maximum ratio ChlNmax 3 mg Chl (mmol N)-1 Ward and Follows, 2016
Table A3.2: Physical Model Parameters
Type Parameter Symbol Value Units Reference
Light
Surface Light PARo varies W m2
Brock 1981
seawater
attenuation Ksw 0.04 m
-1 modified from Li et al., 2010
Chl
attenuation Kchl 0.0249 m
-1 (mg Chl m-3)
-1 Li et al., 2010
Grid Parameter
depth of
each box dz 4 m
157
Table A3.3: Initial Conditions
Parameter Symbol Value Units Reference
Nitrate NO3
varies with
temperature mmol NO3 m-3 Lucas et al., 2011
Ammonia NH4 0 mmol NH3 m-3
Carbon (0-50m) C 10 mmol C m-3
Carbon (50-200m) C 5.00E-05 mmol C m-3
Chlorophyll (0-50m) Chl 2 mg Chl m-3
Chlorophyll (50-100m) Chl 1.00E-06 mg Chl m-3
Total Cells X
varies with
carbon cells m-3
0 5 10 15 20 25
Days until equal biomass with control
0
50
100
150
200
250
Trajectories
Initial Conditions 1
Initial Conditions 2
Figure A3.2: Impact of different initial conditions. Two different sets of initial conditions were compared to
the original initial conditions used in the main text. This sensitivity test was performed on all 295 particles
released at 38 . This histogram shows the days that it took for the biomass of the two sensitivity test
conditions to equal that of the original initial conditions. From this, we determined that initial conditions
largely only impact the first 10 days of the trajectories and so those are excluded from analyses.
158
Abstract (if available)
Abstract
Phytoplankton are vital to the global carbon cycle, contributing half of the world's photosynthesis. This process removes carbon dioxide from the atmosphere, helping regulate Earth's climate. However, large-scale climate models struggle to represent phytoplankton dynamics due to their coarse spatial resolution (~10-150km) which smooths spatial variability. Smoothing spatial variability overlooks the nonlinear response of phytoplankton growth and biomass accumulation rates to environmental changes. This dissertation investigates the critical bio-physical interactions that occur in the ocean that govern the growth of phytoplankton and resulting biomass accumulation. I focus on the timescales ranging from hourly to seasonal dynamics and spatial scales ranging from 1km to the entire Southern Ocean. To accurately determine the bio-physical response to environmental fluctuations, it is essential to adopt the phytoplankton’s reference frame. As phytoplankton are largely passive drifters, they experience environmental variability in the “along-flow” reference frame (Lagrangian reference frame), which is distinctly different from the fixed-point, Eulerian reference frame. Using a combination of Lagrangian based methods and biogeochemical models, I quantify how phytoplankton growth and biomass accumulation respond to varying rates and magnitudes of temperature changes over hours to seasons. I also unravel the complex interactions between advection/diffusion and biomass accumulation at small scales (hourly, ~1km). I show that moderate temperature changes (3-6 oC) within 7-21 days cause lasting decreases in community growth rates, a phenomenon not captured by typical, large-scale model, growth dynamics. Furthermore, resolving fine-scale physics reveals that biomass accumulation is driven by a complex interplay between mixing and non-linear growth dynamics. These findings challenge existing hypotheses on biomass accumulation and highlight the intricate dynamics governing phytoplankton populations. This improved understanding provides a basis for refining large-scale models, ultimately enabling more accurate predictions of the ocean's role in the global carbon cycle and climate regulation.
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
The dynamic regulation of DMSP production in marine phytoplankton
PDF
Future impacts of warming and other global change variables on phytoplankton communities of coastal Antarctica and California
PDF
Thermal diversity within marine phytoplankton communities
PDF
Iron-dependent response mechanisms of the nitrogen-fixing cyanobacterium Crocosphaera to climate change
PDF
Flowstone ideograms: deciphering the climate messages of Asian speleothems
PDF
Concentration and size partitioning of trace metals in surface waters of the global ocean and storm runoff
PDF
The impact of mesoscale and submesoscale physical processes on phytoplankton biomass, community composition, and carbon dynamics in the oligotrophic ocean
PDF
Taking the temperature of the Common Era: statistics, patterns and dynamical insights
PDF
Investigations on marine metal cycling through a global expedition, a wildfire survey, and a viral infection
PDF
Spatial and temporal investigations of protistan grazing impact on microbial communities in marine ecosystems
PDF
New insights into glacial-interglacial carbon cycle: multi-proxy and numerical modeling
PDF
Carbonate dissolution at the seafloor: fluxes and drivers from a novel in situ porewater sampler
PDF
Ecological implications of daily-to-weekly dynamics of marine bacteria, archaea, viruses, and phytoplankton
PDF
Thermal acclimation and adaptation of key phytoplankton groups and interactions with other global change variables
PDF
Abrupt climate variability during the last glacial and the late Holocene from the western tropical Pacific perspective
PDF
The geobiological role of bioturbating ecosystem engineers during key evolutionary intervals in Earth history
PDF
Seeing the future through the lens of the past: fusing paleoclimate observations and models
PDF
Environmental controls on alkalinity generation from mineral dissolution: from the mineral surface to the global ocean
PDF
On the hunt for cryptic thrusts: hot rocks, low geothermal gradients, and elusive structures in the eastern Great Basin, Nevada
PDF
Discerning local and long-range causes of deoxygenation and their impact on the accumulation of trace, reduced compounds
Asset Metadata
Creator
Zaiss, Jessica
(author)
Core Title
The impact of Lagrangian environmental variability on the growth of phytoplankton
School
College of Letters, Arts and Sciences
Degree
Doctor of Philosophy
Degree Program
Geological Sciences
Degree Conferral Date
2024-08
Publication Date
10/02/2024
Defense Date
08/29/2024
Publisher
Los Angeles, California
(original),
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
biophysical interactions,Lagrangian analysis,OAI-PMH Harvest,phytoplankton,submesoscale,variability
Format
theses
(aat)
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Levine, Naomi (
committee chair
), Feakins, Sarah (
committee member
), Hutchins, David (
committee member
), John, Seth (
committee member
)
Creator Email
jzaiss@gmail.com,zaissbow@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-oUC11399BEEE
Unique identifier
UC11399BEEE
Identifier
etd-ZaissJessi-13570.pdf (filename)
Legacy Identifier
etd-ZaissJessi-13570
Document Type
Dissertation
Format
theses (aat)
Rights
Zaiss, Jessica
Internet Media Type
application/pdf
Type
texts
Source
20241002-usctheses-batch-1216
(batch),
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright.
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Repository Email
cisadmin@lib.usc.edu
Tags
biophysical interactions
Lagrangian analysis
phytoplankton
submesoscale
variability