Assimilating bio-optical glider data during a phytoplankton bloom in the southern Ross Sea

18 Pages·2017·5.25 MB·English

Biogeosciences,15,73–90,2018
https://doi.org/10.5194/bg-15-73-2018
©Author(s)2018.Thisworkisdistributedunder
theCreativeCommonsAttribution4.0License.
Assimilating bio-optical glider data during a phytoplankton bloom
in the southern Ross Sea
DanielE.Kaufman1,MarjorieA.M.Friedrichs1,JohnC.P.Hemmings2,a,andWalkerO.SmithJr.1
1VirginiaInstituteofMarineScience,CollegeofWilliam&Mary,GloucesterPoint,VA,USA
2WessexEnvironmentalAssociates,Salisbury,UK
anowat:MetOffice,Exeter,UK
Correspondence:DanielE.Kaufman([email protected])andMarjorieA.M.Friedrichs([email protected])
Received:21June2017–Discussionstarted:30June2017
Revised:22November2017–Accepted:22November2017–Published:4January2018
Abstract.TheRossSeaisaregioncharacterizedbyhighpri- quencyconsistentwithavailablesatellite-deriveddataresults
mary productivity in comparison to other Antarctic coastal in 40% lower carbon export, primarily resulting from opti-
regions,anditsproductivityismarkedbyconsiderablevari- mized rates generating more slowly sinking diatoms. This
ability both spatially (1–50km) and temporally (days to analysishighlightstheneedforthestrategicconsiderationof
weeks). This variability presents a challenge for inferring the impactsof datafrequency, duration, andcoverage when
phytoplankton dynamics from observations that are limited combining observations with biogeochemical modeling in
in time or space, which is often the case due to logistical regionswithstrongmesoscalevariability.
limitationsofsampling.Tobetterunderstandthespatiotem-
poral variability in Ross Sea phytoplankton dynamics and
to determine how restricted sampling may skew dynamical
interpretations, high-resolution bio-optical glider measure- 1 Introduction
mentswereassimilatedintoaone-dimensionalbiogeochem-
ical model adapted for the Ross Sea. The assimilation of Phytoplankton blooms in the Ross Sea are responsible for
datafromtheentireglidertrackusingthemicro-geneticand some of the highest rates of productivity in the South-
local search algorithms in the Marine Model Optimization ern Ocean (Arrigo et al., 2008), and yet the phytoplankton
Testbed improves the model–data fit by ∼50%, generating assemblage exhibits considerable spatiotemporal variabil-
rates of integrated primary production of 104gCm−2yr−1 ity (DiTullio and Smith, 1996; Hales and Takahashi, 2004;
and export at 200m of 27gCm−2yr−1. Assimilating glider Smithetal.,2010).Thisheterogeneity,andthespatialortem-
data from three different latitudinal bands and three differ- porallimitationsofobservationsduetologisticalchallenges
entlongitudinalbandsresultsinminimalchangestothesim- of sampling, may affect the inferred phytoplankton dynam-
ulations, improves the model–data fit with respect to unas- ics and produce biases in productivity or export estimates.
similateddataby∼35%,andconfirmsthatanalyzingthese The magnitude of the underlying ecosystem variability that
glider observations as a time series via a one-dimensional contributes to these potential biases is not well understood,
model is reasonable on these scales. Whereas assimilating nor is it well known how the use of different observational
the full glider data set produces well-constrained simula- platforms in the Ross Sea might affect the inferred dynam-
tions, assimilating subsampled glider data at a frequency ics. Acquiring data with an appropriate resolution is impor-
consistentwithcruise-basedsamplingresultsinawiderange tant for assessing phytoplankton variability in the Ross Sea
ofprimaryproductionandexportestimates.Theseestimates (HalesandTakahashi,2004).
depend strongly on the timing of the assimilated observa- Overthepastseveraldecades,biogeochemistryintheRoss
tions,duetothepresenceofhighmesoscalevariabilityinthis Sea has been observed by ship and satellite, providing data
region.Assimilatingsurfacegliderdatasubsampledatafre- atdifferenttemporalandspatialresolutions.SinceRossSea
phytoplankton became a focus of scientific research in the
PublishedbyCopernicusPublicationsonbehalfoftheEuropeanGeosciencesUnion.
74 D.E.Kaufmanetal.:Assimilatingbio-opticalgliderdata
late1970s,watercolumnmeasurementshaveprimarilycome son et al., 1995, 1996; Matear, 1995; Fennel et al., 2001;
fromresearchvessels(e.g.,El-Sayedetal.,1978;Smithand Friedrichs, 2002; Schartau and Oschlies, 2003; Hemmings
Nelson,1985;Vaillancourtetal.,2003).Typically,sampling etal.,2004;Bagniewskietal.,2011;Doronetal.,2013;Xiao
stationsareseparatedbytensofkilometers(HalesandTaka- and Friedrichs, 2014a, b; Melbourne-Thomas et al., 2015;
hashi, 2004), and although vessels may return to resample Song et al., 2016; Gharamti et al., 2017; Schartau et al.,
astation,theytypicallydonotreturnmorethanonceortwice 2017).Thisentailsasmoothingoroptimizationprocedure,in
inasingleyear.Duringthe1990s,theuseofremotesensing whichelementsofthemodelareadjustedtominimizediffer-
was expanded to look more closely at the Ross Sea bloom encesbetweenthemodeloutputandtheobservations.Typi-
(ArrigoandMcClain,1994),andsatelliteretrievalshavecon- cally,anaggregatemeasureofthedifferencesbetweenobser-
tinuedtoprovidevaluableinsightsintocharacteristicsofthe vationsandmodeloutputisprovidedbycalculationofacost
phytoplankton assemblage (Arrigo et al., 1998; Arrigo and function,definedasthemodel–datamisfit,andanoptimiza-
van Dijken, 2004; Peloquin and Smith, 2007; Schine et al., tionalgorithmsearchesformodelparametersthatminimize
2015).Satelliteobservationsofferasynopticviewofspatial thevalueofthiscostfunction.
regionsatfrequenciesthatarewithinthetimescaleofbiolog- Inthisstudy,dataassimilationisusedtoobtainanoptimal
icalchanges(e.g.,growth);however,thepresenceofseaice representationofRossSealowertrophiclevels.Specifically,
and clouds often obscures remote-sensing measurements in observationsfromanautonomousgliderareassimilatedinto
theRossSea(Arrigoetal.,1998). a biogeochemical model of the Ross Sea (Kaufman et al.,
On the mesoscale (days–weeks, 1–10km), gliders are 2017a)tobetterunderstandthespatialandtemporalvariabil-
a relatively new and effective means to characterize phyto- ityinphytoplanktoninthisregion.Assimilationexperiments
plankton variability, and the development of ice-avoidance also examine how the space and time characteristics of ob-
algorithms has enabled the use of gliders in the Ross Sea servational sampling frequency impact the ability of obser-
forthesepurposes.Forexample,agliderequippedwithbio- vationstoproduceoptimalsystemrepresentations.
optical sensors was directed along a section near 76◦40(cid:48)S
in austral summer 2010–2011 and provided valuable esti-
2 Methods
mates of biomass variability on short timescales (Kaufman
etal.,2014).EstimatesofthePOC:Chlratiofromtheglider
2.1 One-dimensionalbiogeochemicalmodel
optical sensors suggested a transition from a Phaeocystis
antarctica(P.antarctica)toadiatom-dominatedassemblage
Numerical experiments were conducted with the Model of
over several days (Kaufman et al., 2014; Thomalla et al.,
EcosystemDynamics,nutrientUtilisation,Sequestrationand
2017). Moreover, Jones and Smith (2017) used glider ob-
Acidification for the Ross Sea (MEDUSA-RS; Kaufman
servations from austral summer 2012–2013 to distinguish
etal.,2017a),aregionallyadaptedversionofMEDUSA-1.0
three phases of the Ross Sea bloom and identified high-
(Yool et al., 2011). Three phytoplankton groups are repre-
frequency(hours)associationsbetweenwind-drivenmixing
sented in the MEDUSA-RS model: colonial P. antarctica,
and biomass. A perennial challenge when using glider data
solitary P. antarctica, and diatoms. Phytoplankton growth
(as well as ship-based data), however, is separating the ef-
in the model is temperature dependent as well as limited
fectsoftimeandspace(Kaufmanetal.,2014;Little,2016).
by light and nutrient availability. Colonial P. antarctica, di-
Numerical models are another approach for examining
atoms,anddetritusallsinkatdistinctrates.Themodelhan-
phytoplanktonvariabilityintheremoteRossSea,providing
dlesthesinkingoflargedetritalparticlesimplicitlyasafast-
an effective means for coordinating knowledge and under-
sinkinggrouptoavoidissuesrelatedtothescaleofthemodel
standingtheunderlyingsystemcomplexities(Leonelli,2009;
timestepandtoavoidtheneedforanadditionaltracer.Abal-
Vallverdú, 2014). Furthermore, numerical simulations offer
last scheme is used to allow inorganic materials to “pro-
thepossibilityforexperimentalmanipulationsthatwouldbe
tect”avariablefractionofthesinkingorganicmaterialfrom
impractical or impossible in the real system. Such manip-
degradation. The model is configured to focus on dynamics
ulations were implemented in the scenario experiments de-
within the euphotic zone with a vertical resolution of 5m
scribed by Kaufman et al. (2017a) to investigate how pro-
from the ocean surface to 200m. A full description of the
jectedclimatechangesmightalterthedynamicsofthephy-
model and its setup within the Marine Model Optimization
toplanktonassemblage.Theseexperimentsshowedthatear-
Testbed(MarMOT;HemmingsandChallenor,2012),aswell
lieravailabilityoflowlightresultingfromseaicereduction
asthephysicalforcingsderivedfromgliderobservations,are
was the primary driver of projected increases in production
documentedinKaufmanetal.(2017a,b).
andexportandcompositionchangeoverthenextcentury.
Dataassimilation,whichreferstomethodologiesthatsys- 2.2 Dataforassimilation
tematically combine a mathematical model with observa-
tions, is often used in biogeochemical applications (Hof- In situ observations used for the assimilation experiments
mann and Friedrichs, 2001, 2002) to improve estimates of came from an iRobot Seaglider equipped with a Wet
model parameters that are frequently poorly known (Law- Labs ECO Puck sensor and are available in the Biologi-
Biogeosciences,15,73–90,2018 www.biogeosciences.net/15/73/2018/
D.E.Kaufmanetal.:Assimilatingbio-opticalgliderdata 75
strumenterroranderrorofrepresentativeness)andnonpara-
metricerrorsinthesimulationsuchasforcingerrors(Schar-
tauetal.,2017).Theuseofdifferentcostfunctionweighting
schemes in plankton modeling, including the characteristic
scaletechnique,isexploredinmoredetailbyHemmingsand
Challenor(2012).
2.4 Costfunctionminimization
Model parameters were optimized in MarMOT by finding
theminimumofthecostfunction(Sect.2.3)throughacom-
binationofthemicro-geneticalgorithm(µGA)andPowell’s
non-gradientdirectionsetalgorithm.TheµGArunsfirstand
identifiessetsofparametervaluesthatproducelow-costval-
ues; this is achieved by “evolving” a population of various
Figure1.SouthernRossSeashowingtransectlocationswherethe parametersetsoversuccessiveiterations,calledgenerations.
glider was at the surface. The color of each glider dive indicates The low-cost parameter sets identified by the µGA are then
thedate.Bathymetriccontoursareshownat200mintervals,asob- used as starting points for the direction set method, which
tainedfromthebedmap2bathymetricdata(Fretwelletal.,2013). performssuccessivelinearsearchestoidentifynearbylower
costsolutions.
Genetic algorithms, including the µGA, are a subtype of
cal and Chemical Oceanography Data Management Office
computational methods known as evolutionary algorithms,
datarepository(Smith,2005).Gliderdivesfrom22Novem-
so-calledbecauseoftheirinspirationfrom,andsimilarityto,
ber 2012 to 1 February 2013 covered a horizontal area
biologicalevolution.Describedusingthisanalogy,agenetic
spanning76.83–77.44◦Sand168.9–171.97◦E(Fig.1).Data
algorithm procedure modifies a population of candidate so-
spanningtheupper200mofthewatercolumnwerebinned
lutions over successive generations by variation and selec-
by means into hourly, 5m vertical bins. Concentrations
tion processes to converge on a single solution or solution
of chlorophyll (Chl) and particulate organic carbon (POC)
area. Genetic algorithms (GAs) have several advantages for
were derived, respectively, from fluorescence and optical
optimization,includingtheirintrinsicparallelism,suitability
backscattercountsmeasuredbythesensorandconvertedus-
for systems with multiple local minima, and their general-
ingregressionequations(Kaufmanetal.,2017a).Thesebio-
izability (Bajpai and Kumar, 2010; Ward et al., 2010). The
opticalquantitieswereusedforcalculatingmodel–datamis-
µGAusesthreestepstotransitionfromonegenerationtothe
fitsduringassimilation(Kaufmanetal.,2017c).
next,describedfollowingthebiologicalanalogyasselection,
crossover, and resampling (Krishnakumar, 1990; Cˇrepinšek
2.3 Costfunction
et al., 2013). An advantage of the µGA is its reduced risk
The “cost function” (J), defined as a measure of misfit be- ofprematureconvergence,resultingfromre-initializingafter
tweenaparticularmodelsimulationandobservationaldata, each convergence, and generating new random populations
iscomputedasaweightedaverageofthesquareddifferences while maintaining the best fit individual from the previous
betweensimulatedandobservedvalues: set(Schmitt,2001).
In this µGA implementation, optimizations begin with
!
J = 1 XN 1 (cid:0)x −y (cid:1)2+ 1 (cid:0)x −y (cid:1)2 , apopulationoffiveindividualparametersetsrandomlygen-
N σ2 i,chl i,chl σ2 i,poc i,poc erated for the first µGA generation. The constituent pa-
i=1 chl poc
rameter values are selected randomly from within a pre-
whereN isthenumberofobservationpoints,x isthesim- determinedrangeofallowablevalues(Sect.2.5.1).Aneval-
i
ulated value of either chlorophyll or POC at the ith obser- uationofthecostfunctionforeachmodelsolutionindicates
vation point, and y is its observed value; σ is the standard the“fitness”ofeachindividual.Abinarytournamentproce-
i
deviation(SD)ofthespecificobservationsetassimilatedin dure is then followed to select parents from this population
aparticularexperiment.UsingtheSDoftheobservationsto for the next generation. The most-fit individuals (i.e., those
define a characteristic scale of variation for each variable is withthelowestcostfunctionvalues)arepairedwithonean-
a technique used in previous studies (e.g., Friedrichs et al., otherandundergorecombinationofthebitsrepresentingpa-
2006;XiaoandFriedrichs,2014a,b).Itisdesignedtoweight rametervalues.Aftereachgeneration,theproportionofbits
therelativemisfitcontributionofeachvariableappropriately differing from those of the fittest individual is calculated to
whenthereareinsufficientdatatodefineacomprehensiveer- determinewhetherthepopulationcanbedeemedconverged
rormodel.Suchamodelwouldrequirereliableinformation (though this does not necessarily indicate closeness in pa-
abouttheuncertaintyassociatedwithobservationerrors(in- rameterspace).Afterthethresholdforconvergencehasbeen
www.biogeosciences.net/15/73/2018/ Biogeosciences,15,73–90,2018
76 D.E.Kaufmanetal.:Assimilatingbio-opticalgliderdata
achieved,thepopulationisre-initializedtorandomindividu- notexpressedasfractionswerelog-transformedforsampling
als,althoughthefittestindividualismaintained.TheµGAis purposes.
terminateduponthefirstconvergenceoccurringafteramin-
imumnumberofgenerationshasbeenreached. 2.5.2 Sensitivityanalysis
Once convergence has been achieved after a minimum
number of µGA generations, Powell’s non-gradient direc- Parameters to which model outputs are highly sensitive are
tion set algorithm performs a local search using the µGA important and useful to optimize. In contrast, it is futile to
solutions as starting points. The direction set method per- optimize parameters to which model outputs of interest are
formssequentialminimizationsiniterativedirections,updat- not sensitive; no amount of varying these parameters will
ing the search direction after each iteration (Powell, 1964; result in improved model performance. Therefore, the first
Pressetal.,1992).AlthoughtheµGAiswellsuitedforglobal criterion used to designate a parameter as optimizable was
searchproblemspartlybecauseofitsstochasticity,Powell’s thesensitivityofmodeloutputstothevaluesofthatparam-
directionsetalgorithmiswellsuitedtosearchingforalocal eter.Modelsensitivitieswereevaluatedforassimilatedvari-
optimum. Brent’s method, which combines root-bracketing ables (Chl and POC) and carbon fluxes of interest (primary
withsecantandinversequadraticinterpolation(Brent,1973), production (PP) and carbon export at 200m). To quantify
is used to numerically locate cost minima between neigh- the sensitivities of these outputs to each of the 80 parame-
boring function evaluations along each direction identified ters in the model, a series of runs were conducted follow-
by the Powell algorithm. The direction set algorithm stops ing the approach of Hemmings et al. (2015). Each run used
when a cost function minimum is located or when a maxi- auniquesampleofparametervaluesdrawnfromwithinthe
mumnumberofiterationsisreached.Theoptimizedparam- specifiedparameterranges(Sect.2.5.1)usingaLatinhyper-
eter values are those that generated the cost function mini- cube.Thisapproachprovidesmoreevencoverageofthepa-
mum. rameterspace thanMonte Carlosamplingmethods thatcan
resultinclusteredvaluesandunsampledregions(Appendix
2.5 Selectionofparameterstobeoptimized A). One thousand values were drawn from sequential inter-
valsthroughouttherangeforeachparameter.Usingthistech-
Ideally, optimal values are identified for all parameters in nique,uniqueparametersetswereconstructedsuchthatover
a model, however, uncertainty in the parameter estimates the course of all runs, the full range of values for each pa-
from an algorithmic optimization increases as the num- rameterwasrepresented.
ber of parameters included in that optimization increases The model was run 1000 times, each time using one of
(Friedrichsetal.,2007;Wardetal.,2010).Althoughtheop- the unique parameter sets resulting from Latin hypercube
timization of more parameters generally lowers the assimi- sampling of the full parameter space. Sensitivity was quan-
latedcost,theincreasingpotentialforequifinalitywithmore tifiedbyevaluatingtheamountofvarianceintheoutputdi-
parametersmeanstheoptimizationmayfindequivalentlow- agnostics explained by each parameter (i.e., by computing
cost solutions with substantially different parameter values. thecoefficientofdetermination(r2)betweeneachparameter
Therefore, before assimilating observations and optimizing andeachofthefouroutputvariablesofinterest;Fig.2).All
parameters, a subset of “free” or “optimizable” model pa- four model outputs (Chl, POC, PP, and export) were most
rametersmustbechosen.Inthisstudy,theparameterstobe sensitive (r2≥0.01) to attenuation of blue-green light by
optimizedareselectedthroughathree-stepprocess:defining phytoplanktonpigments,diatommaximumgrowthrate,and
arangeofpermittedvaluesforeveryparameter(Sect.2.5.1), C:Chl ratio for solitary P. antarctica. Three additional pa-
identifyingtheparameterstowhichmodeloutputsaremost rameters (maximum growth rate of P. antarctica colonies,
sensitive(Sect.2.5.2),andevaluatinghowmanyofthesesen- maximum growth rate of solitary P. antarctica, and micro-
sitiveparameterscanbereasonablyoptimizedwhenassimi- zooplanktonmaximumgrazingrate)exhibitedr2≥0.01for
latingtheavailabledata(Sect.2.5.3).Initialvaluesforeach bothchlorophyllandPOC.The21parameterswithr2≥0.01
parameter,priortotheassimilation,weresettovaluesiden- (Fig.2)wereselectedforfurtherevaluation(Sect.2.5.3).
tifiedinKaufmanetal.(2017a).
2.5.3 Usingtwinexperimentstoselectoptimizable
2.5.1 Parameterranges subset
Upper and lower bounds of the allowable range for each After selecting the 21 potentially optimizable parameters,
free parameter were defined loosely following Hemmings numericaltwinexperiments(NTEs)wereconductedtoiden-
etal.(2015).Boundsweresettobegeometricallysymmetric tifyanoptimizablesubsetbyevaluatingtheextenttowhich
(factor of 4 for rates; factor of 5 for half-saturation concen- known values of sensitive parameters could be recovered
trations) around the initial values. For fractional parameter given the data available for assimilation. The implementa-
values,limitsweresetto±0.25theirinitialvalues,although tion of NTEs involves four primary steps (Hofmann and
not allowed to exceed 0.05 or 0.95. Ranges for parameters Friedrichs, 2001). First, the chosen model is run forward in
Biogeosciences,15,73–90,2018 www.biogeosciences.net/15/73/2018/
D.E.Kaufmanetal.:Assimilatingbio-opticalgliderdata 77
Figure2.VarianceexplainedinmodeloutputsbyparametersduringsensitivitytestsusingLatinhypercubesamplingofparameterspace.
Onlyparameterswithatleastoner2valuegreaterthanorequalto0.01(verticaldottedline)areshown.
time to create a simulation using a known, “true” parame- tions(Wardetal.,2010).Infact,ratherthanbeingafunction
terset.Second,outputfromthissimulationissubsampledto oftheoptimizationalgorithm,itisdependentontheavailable
create a so-called “synthetic” data set. Third, the synthetic dataandthedesignofthecostfunction(LöptienandDietze,
dataset is then assimilated to optimize model parameters. 2015).Alargerorricherobservationsetcanhelptoconstrain
Fourth, the optimized parameter set is compared to the true moreparameters.Theimpactofcostfunctiondesignismore
parameterset.Theassimilationissuccessfuliftheoptimized complicated because an improved cost function may allow
values recover the true parameters used to generate the as- for greater uncertainty in the observations and/or nonpara-
similatedsyntheticdata. metric uncertainty in the simulation, leading to weaker but
There is a limit to the number of parameters that can be morerealisticconstraintsontheparameters(Hemmingsand
independently constrained by the available observations be- Challenor,2012).
causevaryingdifferentparameterscanoftenhavesimilaref- Theprocedurefollowedherefordeterminingthesubsetof
fects on the cost function. Optimizing a larger set increases optimizable parameters is similar to that used by Friedrichs
the potential for correlation between the effects of different etal.(2007).First,areferencesimulationwasgeneratedus-
parameters,reducingthealgorithm’seffectivenessinidenti- ing the initial parameter set, and chlorophyll and POC es-
fying unique optimal parameter sets. This, combined with timates from this reference simulation were subsampled to
the increased potential for over-fitting associated with the generateasyntheticdataset.Startingwithaparameterspace
greater model degrees of freedom, can reduce the ability of defined by the set of 21 parameters deemed sensitive in the
anoptimizedmodeltoreproduceindependentdata(Matear, Latin hypercube tests (Fig. 2), a series of sequential NTEs
1995; Friedrichs et al., 2007; Xiao and Friedrichs et al., wasthenperformedwithaprogressivelysmallernumberof
2014b).Thelimitationonthenumberofoptimizableparam- optimizedparameters:aftereachNTE,theoptimizedparam-
eters applies to both µGA and variational adjoint optimiza- eter that was most different from its “true” value was re-
www.biogeosciences.net/15/73/2018/ Biogeosciences,15,73–90,2018
78 D.E.Kaufmanetal.:Assimilatingbio-opticalgliderdata
only minor differences between assimilating data from the
upper50mvs.thefullupper200m.)Observationsfromdif-
ferentspatialareasoftheglidertrackwerealsoassimilated.
Observations from the glider track were divided into three
latitudinal bands (northern, central, and southern bands) as
wellasintothreelongitudinalareasconstitutingeastern,cen-
tral,andwesternbands.Gliderdatafromeachofthesethree
latitudinalandlongitudinalbandswereassimilatedinExper-
iments 1band1c,respectively(Table2,Fig.4),resultingin
threecostfunctionsforeachoftheseexperiments.
2.6.2 Experiment2
The second set of experiments investigated the assimila-
tionofdataatdifferentresolutionsmimickingdifferentdata
sources. In Experiment 2a, glider data were subsampled at
∼12h intervals (Table 2). The subsampling was repeated
12times,witheachiterationoffsetfromthepreviousby+1,
2,3,...11h,togenerateaseriesof12gliderobservationsets.
The assimilation of these 12 time series yields the “Glider
Figure 3. (a) Minimum costs and (b) normalized parameter val-
Assimilation”case.InExperiment2b,gliderdataweresub-
uesinnumericaltwinexperiments,illustratingthattheassimilation
sampled at a reduced temporal resolution similar to cruise
procedureisunabletosuccessfullyrecoverthetrueparameterval-
sampling (Table 2). Sampling during cruise missions often
ueswhenmorethaneightparametersareoptimized.Onedatapoint
takesplaceforafewdaysinonelocationbeforemovingelse-
inthreeoftheexperiments(19,20,and21)exceedstheyaxisupper
limitinpanel(b). where, and the ship sometimes returns to the first location
after a number of weeks. To roughly mimic this sampling
pattern,dailyverticalprofiles(againdownto50m)wereas-
movedfromtheoptimizableparameterset.Thus,aftereach similated for 3 days in a row, starting from the first day of
NTE the number of optimized parameters was reduced by availablegliderdata(22November),andthen3daysofdata
one.TheseriesofNTEswasevaluatedtoidentifythelargest wereassimilated2weekslater.Shiftingthispatternforward
parameter set for which the original parameter values were 1weekatatimegeneratedaseriesofeightcruise-basedob-
recoverableandforwhichthecostfunctionremainedessen- servationsetsforassimilationinthis“Cruise-BasedAssimi-
tially 0. From this analysis (Fig. 3), it was determined that lation”case.InExperiment2c,gliderdatawereassimilated
optimizingeightparameterswouldbeideal(Table1)because only from the upper 5m surface layer to produce a data set
valuesoftheseeightparameterswererecoveredmuchbetter resemblingsatellite-deriveddata.Thesedatawerethensub-
than larger parameter sets and model–data misfit (cost) re- sampledat2-weekintervals,torepresenttypicaldatareturn
mainedlow. fromremote-sensingobservationsofoceancolorintheRoss
Sea,wheretheavailabilityofsatelliteimageretrievalisfre-
2.6 Assimilationexperiments quently limited by excessive, though variable, cloud cover
(ArrigoandvanDijken,2004).The2-weeksubsamplingpat-
The µGA optimization procedure was used to assimilate terncoveredtheentireperiodofgliderdata(22November–
glider data in two sets of experiments that explored aspects 1 February) and was sequentially shifted forward 1 day at
ofspatiotemporalvariabilityanddataavailability.Estimates a time to generate a series of 14 satellite-based observation
ofdepth-andtime-integratedPPandtime-integratedcarbon sets for assimilation in this “Satellite-Based Assimilation”
exportat200mwerecomputedfromthefullmodelsimula- case(Table2).
tionineachexperiment.
2.7 Predictivecostassessment
2.6.1 Experiment1
In addition to the assimilative cost (J ) calculated during
A
The first set of experiments examined the differences in theoptimizationprocedureusingassimilateddata,apredic-
modelsimulationsresultingfromassimilatingChlandPOC tive cost (J ) was calculated to assess model–data misfit
P
datafromdifferentspatialregions.InExperiment1a,glider computed using the unassimilated data in each experiment.
observations were assimilated from the upper 50m of the Because predictive costs represent model–data misfit from
fulltemporalandspatialdomain,referredtohereafterasthe unassimilateddataonly(Friedrichsetal.,2006;Wardetal.,
“Full Assimilation” case (Table 2). (Comparisons showed 2010),itisanobjectivemeasureoftheskillofanoptimized
Biogeosciences,15,73–90,2018 www.biogeosciences.net/15/73/2018/
D.E.Kaufmanetal.:Assimilatingbio-opticalgliderdata 79
Table1.Eightparametersoptimizedinthisanalysis.
Parametername Initialvalue Bounds
(Kaufmanetal.,2017a) (lower,upper)
Diatommax.growthrateat0◦C 0.375(d−1) 0.09375,1.5
P.antarcticasolitarycellsC:Chlratio 30(gCgChl−1) 7.5,120
P.antarcticacoloniesmax.growthrateat0◦C 0.5(d−1) 0.125,2
P.antarcticasolitarycellsmax.growthrateat0◦C 0.5(d−1) 0.125,2
DiatomC:Chlratio 150(gCgChl−1) 37.5,600
Fastdetritussinkingfractionofdiatomlosses 0.75 0.5,0.95
P.antarcticacoloniesmax.sinkingrate 20(md−1) 5,80
P.antarcticacoloniesC:Chlratio 40(gCgChl−1) 10,160
Figure4.Locationsofgliderobservationsassimilatedin(a)Experiment1b–latitudinalbands–and(b)Experiment1c–longitudinalbands.
Colorsrepresentthethreespatialbandsofdataassimilated.
modelinreproducingobservationsatdifferentpointsintime isloweranddiatomcarbonishigherinDecemberandearly
orspace(Greggetal.,2009).Inthiscase,theaimoftheseex- January (Fig. 6b and d). Compared to the No Assimilation
perimentsistoassesstheskillofeachoptimizedsimulation case, PP is only slightly lower (7%), whereas export flux
regardless of which subset of the available data is assimi- is nearly 50% higher (Table 3; Fig. 7). Compared to their
lated.Bycomputingthemeanandmedianpredictivecostfor initial values, colonial P. antarctica parameters change the
eachexperiment(otherthantheFullAssimilationcase),the mostasaresultoftheoptimization,withreductionsbetween
skill of the resulting simulations can be compared directly 40and70%forthecolonialP.antarcticamaximumgrowth
withoneanother. rate, maximum sinking rate, and C:Chl ratio (Table 4). In
contrast,thediatommaximumgrowthrateandC:Chlratio
increased(∼10and20%,respectively).
3 Results Chlorophyll and POC time series exhibit only minor dif-
ferences between latitudinal band experiments when data
3.1 Experiment1 from the northern, central, and southern sections are as-
similated independently (Fig. 5c and d) or when data from
Theassimilationofthegliderdataoverthefulltemporaland
the eastern, central, and western sections are assimilated
spatialdomain(FullAssimilationcase)improvesthemodel–
(Fig.5eandf).Specifically,theoptimalsimulationsforChl
data fit of both Chl and POC (Fig. 5a and b) and reduces
andPOCexhibitsimilarseasonalcyclesacrossthethreelat-
thecostbynearlyhalf(47%)comparedtotheapriorisimu-
itudinal and longitudinal bands, with only slightly higher
lationwithoutassimilation(NoAssimilationcase;Table3).
Chl and POC concentrations when assimilating data from
Average Chl and POC concentrations in the upper 50m are
the southern band (Fig. 5c and d) and higher Chl from the
both slightly lower (8 and 12%, respectively) in the opti-
westernband(Fig.5eandf).Meancostsaremuchlowerfor
mized simulation. The contribution of each phytoplankton
the latitudinal and longitudinal experiments than for the No
grouptototalchlorophyllremainssimilartotheNoAssimi-
Assimilationcase,andonlyslightlyhigherthantheFullAs-
lationcase(Fig.6aandc),butcolonialP.antarcticacarbon
www.biogeosciences.net/15/73/2018/ Biogeosciences,15,73–90,2018
80 D.E.Kaufmanetal.:Assimilatingbio-opticalgliderdata
Table2.Spatiotemporalresolutionofglider-basedobservationsofChlandPOCassimilatedforeachexperiment.
Experiment Depth Temporalresolution Spatialarea(s)
(m)
Expt1a:FullAssimilation 0–50 Hourly Fullglidertrack
Expt1b:LatitudinalAssimilation 0–50 Hourly North,central,southlatitudinalbands
Expt1c:LongitudinalAssimilation 0–50 Hourly East,central,westlongitudinalbands
Expt2a:GliderAssimilation 0–50 ∼twiceperday,separatedata Fullglidertrack
minimumof12h.
Expt2b:Cruise-BasedAssimilation 0–50 3daysinarow,andthenanother3 Fullglidertrack
consecutivedays2weekslater
Expt2c:Satellite-BasedAssimilation 0–5 1dayevery2weeks Fullglidertrack
Figure5.Upper50mmeanconcentrationsof(a,c,e)Chland(b,d,f)POCforvariousexperimentsassimilatingthefullgliderandfrom
differentspatialareas(Table2):(a,b)Experiment1a,(c,d)Experiment1b–latitudebands–and(e,f)Experiment1c–longitudebands.For
reference,gliderdata(blacklines)withshading(gray)representing1SD(fromtheupper50m)areincludedineachpanel.Coloredboxesat
thetopofeachpanelindicatetimesofassimilatedobservations.
similation case (Table 3). This indicates that data sampled 3.2 Experiment2
from within only one spatial band improved the match be-
tween modeled and observed variables in the unassimilated The assimilation of data subsampled at a frequency 1/12
areasaswell.AverageestimatesofPPandexportinboththe thatoftheoriginalgliderdata(Experiment2a)resultsin12
latitudinalandlongitudinalexperimentsareonlyslightlyless model simulations, all of which are similar to the Full As-
(<5%)thantheFullAssimilationestimate(Fig.7,Table3). similation case, with Chl and POC time series closely fol-
Biogeosciences,15,73–90,2018 www.biogeosciences.net/15/73/2018/
D.E.Kaufmanetal.:Assimilatingbio-opticalgliderdata 81
Figure6.Upper50mmeanconcentrationsofthethreephytoplanktongroupsintermsof(a,c)Chland(b,d)POCfortheNoAssimilation
case(a,b)andtheFullAssimilationcase(c,d).Thegliderdataareshown(blackline)withshading(gray)thatrepresents1SDdaily.
Table3.Depth-andtime-integrated(overthelengthofthesimulation,representingyearlyrates)primaryproduction(PP),carbonexport
fluxat200m,andcostsfortheNoAssimilationrun(cost:0.77),Experiment1,andExperiment2.Costsprovideameasureofthemisfit
betweenaparticularmodelsimulationandobservations,andthecostsshownrepresentmean±1SDofassimilativeruns.Theassimilative
andpredictivecostsarecomputedfromtheassimilatedandunassimilateddata,respectively.
Simulationname PP Export Predictive Assimilation
(gCm−2yr−1) (gCm−2yr−1) cost(JP) cost(JA)
NoAssimilation 111.7 18.8 – –
Expt1a:FullAssimilation 104.2 27.2 – 0.41
Expt1b:LatitudinalAssimilation 101.8±3.3 26.1±2.1 0.49±0.13 0.43±0.14
Expt1c:LongitudinalAssimilation 103.2±2.1 26.9±2.1 0.50±0.10 0.46±0.13
Expt2a:GliderAssimilation 103.7±1.8 27.0±1.2 0.43±0.01 0.43±0.03
Expt2b:Cruise-BasedAssimilation 113.1±22.3 24.8±6.6 1.24±0.95 0.52±0.19
Expt2c:Satellite-BasedAssimilation 114.1±10.7 16.7±2.7 1.04±0.36 0.26±0.16
lowing the observed seasonal pattern (Fig. 8a and b). Mean range from mid-November to early January. The mean pre-
assimilative and predictive costs in the Glider Assimilation dictive cost from this experiment (1.24) is roughly 3 times
case are close to the cost of the Full Assimilation case (Ta- theassimilativecostfortheFullAssimilationcase(0.41)and
ble 3). Mean PP and export estimates are also close to es- 3 times the predictive cost for the Glider Assimilation case
timates from the Full Assimilation case. The mean optimal (0.43;Table3).ThePPestimatesfromtheCruise-BasedAs-
parametervaluesobtainedfromtheGliderAssimilationcase similationcasespanabroadrange(92to156gCm−2yr−1)
aregenerallywithin1SDoftheoptimalvaluesfromtheFull around the Full Assimilation estimate but are generally
Assimilationcase(Table4). higher(Fig.7a).Thisexperimentsimilarlyyieldsaverylarge
The assimilation of data subsampled with a frequency range of export estimates (11 to 33gCm−2yr−1) encom-
typical of cruise observations (Experiment 2b) results in passingtheresultsfromExperiment1(Fig.7b).Optimalpa-
a wide range of solutions, with several Chl and POC time rametervaluesobtainedfromtheCruise-BasedAssimilation
series exhibiting markedly different peak bloom timings casearegenerallylesswellconstrained(higherSDs)thanthe
(Fig. 8c and d). Two of the solutions yield substantially GliderAssimilationcase(Table4).
higher concentrations of POC in November, and Chl peaks
www.biogeosciences.net/15/73/2018/ Biogeosciences,15,73–90,2018
82 D.E.Kaufmanetal.:Assimilatingbio-opticalgliderdata
Table4.Initialparametervalues(NoAssimilation)andoptimalparametervaluesafterconductingtheFullAssimilation,GliderAssimilation,
Cruise-BasedAssimilation,andSatellite-BasedAssimilationexperiments.
Parametername Initial Expt1a Expt2a Expt2b Expt2c
∗
value Full Glider Cruise- Satellite-
∗ ∗
Assimilation Based Based
Diatommax.growthrateat0◦C 0.375 0.40 0.43±0.01 0.42±0.15 0.41±0.09
(d−1)
P.antarcticasolitarycells 30 29.7 25.84±5.16 37.3±26.7 51.5±26.8
C:Chlratio(gCgChl−1)
P.antarcticacolonies 0.5 0.29 0.22±0.10 0.45±0.58 0.29±0.17
max.growthrateat0◦C(d−1)
P.antarcticasolitarycellsmax. 0.5 0.39 0.45±0.06 0.75±0.70 0.79±0.51
growthrateat0◦C(d−1)
DiatomC:Chlratio 150 176.4 166.6±50.17 252.4±164.28 374.86±187.82
(gCgChl−1)
Fastdetritussinkingfractionof 0.75 0.87 0.86±0.05 0.86±0.11 0.62±0.14
diatomlosses
P.antarcticacoloniesmax. 20 10.7 10.1±3.66 20.1±20.5 12.8±9.27
sinkingrate(md−1)
P.antarcticacolonies 40 14.0 14.2±2.29 42.7±41.6 34.3±26.5
C:Chlratio(gCgChl−1)
∗Mean±1SDofassimilativeruns.
Theassimilationofdatasubsampledassatellite-basedob- ineithertheFullAssimilationorCruise-BasedAssimilation
servations from the surface layer (Experiment 2c) results in cases (Table 4). Standard deviations of optimal parameters
Chl and POC concentrations generally higher than the Full fortheSatellite-BasedAssimilationcaseweregenerallysim-
Assimilation case (Fig. 8e and f). The predictive costs are ilartoorlowerthanthosefortheCruise-BasedAssimilation
similaronaveragetothoseoftheCruise-BasedAssimilation case,exceptfortheC:Chlratiofordiatoms,whichproduced
experiment; however, there is less variation (Table 3). The a very high optimal value and was particularly poorly con-
median integrated production is higher (9%) than the Full strained(375±187gCgChl−1;Table4).
AssimilationestimateandtheCruise-BasedAssimilationes-
timate(Fig.7a;Table3);however,therangeofPPestimates
for this Satellite-Based Assimilation case is smaller than 4 Discussion
thosefortheCruise-BasedAssimilationcase(Fig.7a).Most
notably, despite generally higher PP and higher POC con- 4.1 RossSeasimulationresultingfromtheassimilation
centrations, carbon export from the Satellite-Based Assimi- ofgliderdata
lationcaseissubstantiallylower(41%)thantheFullAssim-
Data assimilation is a valuable tool for efficiently utilizing
ilation estimate (Fig. 7b; Table 3). In fact, export estimates
fromindividualrunsinthisexperimentarealllower(−19to limited observational data in remote regions like the Ross
−56%)thantheFullAssimilationestimate(Fig.7b).Again, Sea.Inthisstudy,gliderdataconsistingofbothfluorescence-
derived chlorophyll and backscatter-derived POC were as-
the range of export estimates is smaller for the Satellite-
similated into a one-dimensional marine biogeochemical
BasedAssimilationthanfortheCruise-BasedAssimilation.
model developed for the Ross Sea. Eight ecosystem param-
When assimilating data at a resolution similar to that of
eters,includingphytoplanktongrowthandsinkingratesand
satellite-based observations, mean optimal parameter val-
C:Chlratios,wereoptimizedresultinginasimulationwith
ues were similar to those obtained in the Glider Assimila-
a50%reducedmodel–datamisfit.ThisFullAssimilationrun
tionandCruise-BasedAssimilationcases,withtheexception
yieldedlowerP.antarcticacarbonconcentrationsandhigher
of the fast detritus sinking fraction for diatoms, which was
diatomcarbonconcentrations,resultinginhighercarbonex-
significantly lower in the Satellite-Based Assimilation case
(0.62±0.14) than in the other experiments (Glider Assimi- port compared to those generated by the No Assimilation
lation Case: 0.86±0.05). In contrast to this sinking param- run(Kaufmanetal.,2017a),despiteslightlylowerestimates
of overall annual primary production. Changes in chloro-
eter for mortality from diatoms, the mean maximum sink-
ing rate of colonial P. antarctica in the Satellite-Based As- phyll concentrations of diatoms and P. antarctica were mi-
nor. This Full Assimilation simulation was obtained largely
similation case was not significantly different than its value
via changes in the C:Chl ratios: the colonial P. antarc-
Biogeosciences,15,73–90,2018 www.biogeosciences.net/15/73/2018/