WAP-1D-VAR v1.0: Development and Evaluation of a One-Dimensional Variational Data Assimilation Model for the Marine Ecosystem Along the West Antarctic Peninsula

The West Antarctic Peninsula (WAP) is a rapidly warming region, with substantial ecological and biogeochemical responses to the observed change and variability for the past decades, revealed by multi-decadal observations from the Palmer Antarctica 15 Long-Term Ecological Research (LTER) program. The wealth of these long-term observations provides an important resource for ecosystem modelling, but there has been a lack of focus on the development of numerical models that simulate time-evolving plankton dynamics over the Austral growth season along the coastal WAP. Here we introduce a one-dimensional, data assimilation planktonic ecosystem model (i.e., the WAP-1D-VAR model v1.0) equipped with a variational adjoint and model parameter optimization scheme. We first demonstrate the modified and newly added model schemes to the pre-existing 20 food-web and biogeochemical components of the other ecosystem models that WAP-1D-VAR model was adapted from, including diagnostic sea-ice forcing and trophic interactions specific to the WAP region. We then present the results from model experiments where we assimilate eleven different data types from an example Palmer LTER growth season (October 2002 - March 2003) directly related to corresponding model state variables and flows between these variables. The iterative, data assimilation procedure reduces by 58% the misfits between observations and model results, compared

The wealth of Palmer LTER time-series observations provides an important resource for ecological and biogeochemical modelling, and different types of modelling approaches have been developed to explore the WAP responses to climate change and variability.For instance, an inverse modelling study estimated the steady-state dynamics of the WAP food-web by deriving snapshots of flows among different plankton functional types and higher trophic levels (Sailley et al., 2013).However, there has been a less focus on numerical ecosystem models that simulate time-evolving plankton dynamics over the full Austral growth season along the coastal WAP.Numerical ecosystem models provide estimates of key rate processes for which observations have been less frequently or seldom made compared to frequently measured stocks and rates.Despite its strengths, constructing an ecosystem model is a challenge due to the lack of a priori knowledge on model parameter values and incomplete understanding of ecological processes that should be explicitly presented in the model structure (Ducklow et al., 2008;Murphy et al., 2012).Owing to many observational studies, a more robust, yet still incomplete, databased picture is emerging of WAP food-web interactions and ecosystem dynamics, which could guide a development of the WAP-specific numerical ecosystem model.
Here we introduce a one-dimensional (1-D) variational data assimilation model specific to the coastal WAP (i.e., the WAP-1D-VAR model v1.0) that we develop by adapting an existing biogeochemical-planktonic model of different ocean basins (Friedrichs, 2001;Friedrichs et al., 2006Friedrichs et al., , 2007;;Luo et al., 2010).The WAP-1D-VAR model is compared against the roughly semi-weekly, bio-physical observations over the Austral growth season near Palmer Station on Anvers Island, Antarctica (64.77°S, 64.05°W).The field data record the seasonal variations in the initiation, peak, and termination of phytoplankton blooms and other biogeochemical processes modulated by variations in surface light, mixed layer depth, and sea-ice cover.In the present study, we 1) describe the structure and schemes of the WAP-1D-VAR in great detail, 2) evaluate the model performance and robustness using a variety of quantitative metrics, and 3) discuss the model applicability with regard to capturing the key WAP ecological and biogeochemical features using the data from an example growth season.

Model state variables
The WAP-1D-VAR model v1.0 (Figure 1) is originally derived and modified from data-assimilative, ocean regional test-bed models of the Arabian Sea, the Equatorial Pacific, and the Hawaii Ocean Time-series Station ALOHA (Friedrichs, 2001;Friedrichs et al., 2006Friedrichs et al., , 2007;;Luo et al., 2010).The WAP-1D-VAR model simulates stocks and flows of C, N, and P through 11 different model prognostic state variables.The two size-fractionated phytoplankton compartments (i.e., diatoms and cryptophytes) and the two different zooplankton compartments (i.e., microzooplankton and krill) are separately simulated following the plankton functional types as in Sailley et al (2013) and the observations of the phytoplankton community structure along the coastal WAP.Typically, the coastal WAP is associated with large phytoplankton accumulations dominated by large (> 20 μm) diatoms, but nanoflagellates (< 20 μm) or cryptophytes are also an important component of the food web (Schofield et al., 2017).Mixed flagellates, prasinophytes, and type-4 haptophytes are also found in the region, but we choose to model only diatoms and cryptophytes, in order to avoid too many free (optimizable) parameters associated with each phytoplankton group.The third most dominant species is mixed flagellate but little is known about this taxa in the region and this taxa generally exhibit low interannual variability (Schofield et al. 2017).Functional grazing relationships are defined in which diatoms are consumed by both krill (Euphausia superba) and microzooplankton (mostly ciliates and other protozoa), cryptophytes are consumed by microzooplankton, and microzooplankton are grazed by krill.Other abundant zooplankton taxa in the WAP, such as salps, pteropods, and copepods (Steinberg et al., 2015), are not explicitly simulated in the WAP-1D-VAR model, in part to limit the model complexity and in part because of the limited data constraints on these groups, especially feeding.Higher trophic levels are implicitly represented to close the model.The WAP-1D-VAR model allows for the partitioning between labile dissolved organic matter (LDOM) and semi-labile DOM (SDOM) such that the entire LDOM pool is available but only a limited portion of the SDOM is available for bacterial utilization to account for lower lability of SDOM.Refractory DOM (RDOM) is not explicitly modelled due to its much longer turnover time than labile and semi-labile pools, but some mass flows are included to RDOM from other prognostic model compartments, such as bacteria, krill, and SDOM, to account for loss terms for those state variables.Detritus represents an average particulate organic matter (POM) pool after removing living phytoplankton and bacterial biomasses, and sinking of the detritus pool contributes to particle export flux.The WAP-1D-VAR model explicitly simulates NO3, NH4, and PO4 for inorganic (macro)nutrient compartments, but there is not a separate Fe model compartment or Fe uptake processes, given that even during the peak of the blooms iron is still measurable and iron limitation is absent or occurs only minimally and seasonally in the nearshore Palmer Station area (Carvalho et al., 2016;Sherrell et al., 2018).

Model equations
Here we demonstrate key model processes that are either based on the existing schemes or built as new schemes for the coastal WAP region.The original model schemes are detailed in Supplementary Material of Luo et al (2010).The WAP-1D-VAR model v1.0 simulates biological-physical model processes for a 1-D vertical water column, solving numerically for a discretized version of the time-rate of change for each model state variable.For a generic tracer variable C the time-rate of change equation takes the form (Glover et al., 2011): where z is the depth, w is the vertical velocity (the sum of water motion and gravitational particle sinking), Kz is the turbulent eddy diffusivity, and JC is the biological and biogeochemical net source and sink term for C (Appendix A Equations; Eq.A.2.41-44, A.3.37-40, A.4.53-55, A.5.24-26, A.6.27-29, A.7.4-6, A.8.4-9, A.9.2-4).The physical advection and mixing terms are discussed below in section 2.3 and applied sequentially following the computation of the biological and biogeochemical terms JC using a constant time step of 1 hour.The contributions of the source sink terms JC to the full time rate of change equations are constructed as a series of coupled ordinary differential equations, detailed in Appendix A (sections A1-9), and solved using a second-order Runge-Kutta numerical integration scheme.The WAP-1D-VAR model simulates the dynamics of C, N, and P, but here we only focus on the presentation of the model C dynamics.The cellular molar (e.g., N/C, P/C) quota parameter values of most state variables are fixed (Table 1) and not submitted to the optimization and data assimilation procedure.To first order, most model physiological processes are affected by water temperature, including the maximum growth rates of phytoplankton, bacteria, and zooplankton and basal respiration rates of bacteria and zooplankton.The Arrhenius function is implemented to change these physiological rates as a function of water temperature (Eq.A.1.1).
The net change of phytoplankton (both diatoms and cryptophytes) C biomass is driven by gross growth, DOC excretion, particulate organic carbon (POC) production via aggregation, respiration, and grazing (Eq.A.2.41, A.3.37), the net change of their N and P biomass by gross growth, dissolved organic nitrogen (phosphorus) excretion, particulate organic nitrogen (phosphorus) production, and grazing (Eq.A.2.42-43, A.3.38-39).The net change of their chlorophyll a (Chl) by gross growth, DOM excretion, and grazing (Eq.A.2.44, A.3.40).The WAP-1D-VAR model adapts a phytoplankton growth scheme with flexible stoichiometry, in which phytoplankton cells are allowed to accumulate and store more nutrients under light stress (Bertilsson et al., 2003;Droop, 1974Droop, , 1983;;McCarthy, 1980).The phytoplankton C growth rate is limited by their cellular nutrient quota (Eq.A.2.1-2, A.3.1-2).Modified from Geider et al. (1998), phytoplankton nitrogen uptake decreases when their cellular N/C quota is higher than their reference (Redfield) ratio, but not limited when lower than their reference ratio (Eq.A.2.4, A.2.8, A.3.4, A.3.8).The nitrogen consumption completely ceases when the phytoplankton cellular quota reaches their maximum allowable ratios and is additionally limited by the ambient NO3 and NH4 concentrations with a Monod function (Eq.A.2.10-11, A.3.10-11).NH4 inhibition on NO3 uptake of phytoplankton is modelled by assigning lower k NH4 compared to k NO3 (Table 1).The inhibition term does not exist for PO4.The uptake scheme is similar for PO4 (Eq.A.2.13, A.3.13), but PO4 can be consumed in great excess of current needs (Armstrong 2006).Such luxury uptake is modelled by assigning smaller maximum and minimum P quota, which acts to alleviate P limitation.The maximum photosynthesis rate decreases when the phytoplankton cellular quota is lower than their reference ratio, and approaches zero near their minimum ratio (Eq.A.2.6, A.3.6).The Chl production decreases with lowering photosynthetic active radiation (PAR) and completely ceases in dark (Eq.A.2.14, A.3.14).Phytoplankton release LDOM via passive diffusion of the low molecular weights DOM (e.g., neutral sugars and dissolved free amino acids) with the same cellular elemental ratio as that of phytoplankton (Fogg 1966, Bjørnsen 1988, Biddanda & Benner 1997;Eq. A.2.16-18, A.3.16-18).Phytoplankton also release L-and SDOM actively, in the form of carbohydrate, as 75% of the labile (Eq.A.2.19, A.2.23, A.3.19, A.3.23) and 25% of the semi-labile pools (Eq.A.2.20, A.2.26, A.3.20, A.3.26).This active DOM production enables phytoplankton to adjust their stoichiometry to approach their reference ratio.If cellular organic C is in excess, DOC is released on a time scale of 2 days, and if excess nitrogen (phosphorus), DON (DOP) is released on a time scale of 8 days (Eq.A.2.21-22, A.3.21-22).Diatoms are grazed by both microzooplankton and krill (Eq.A.2.33-40), while cryptophytes are only grazed by microzooplankton (Eq.A.3.33-36).Microzooplankton grazing functions are altered by assigning grazing limitation terms (ϵ) to provide a limit on diatom grazing and route more cryptophytes to microzooplankton (Eq.A.2.33, A.3.33), based on initial modelling attempts where elevated diatom Chl was not simulated due to their much stronger removal by microzooplankton than cryptophytes.In principle, optimization should be able to capture the elevated diatom Chl by adjusting free parameters unless: 1) the right parameters are not adjusted and/or the baseline (non-optimized) parameters need significant adjusting, and/or 2) the model equations are not adequate even with the optimized parameters.In our initial modelling attempts, the model failed to simulate the elevated diatom Chl with varying sets of the model initial parameter values assigned to decouple diatoms from their grazers.Thus, grazing limitation terms (ϵ) are instead assigned to limit microzooplankton grazing on diatoms for modelling purposes, the implementation of which is not strictly based on the ecological evidence of prey switching, or of zooplankton mortality thresholds.
The net change of bacterial biomass is driven by their gross growth respiration (Eq. A.4.25),, and mortality due to viral attack .The WAP-1D-VAR model allows both L-and SDOM as the substrate sources for bacteria, and bacterial nutrient quota lets the lability of SDOM variable for their selective utilization.All the LDOM pool is available, while only a limited portion of the SDOM pool is allowed for bacterial utilization, the degree of which is controlled by an optimizable parameter controlling the relative utilization of SDOM to LDOM,or SDOM lability (i.e.,rSDOC,Eq. A.4.11, Table 1).Bacterial C growth is determined by their cellular quota and available L-and SDOC concentration (Eq.A.4.12-13), in which the growth would be limited if bacterial cellular nitrogen (phosphorus) quota is smaller than their reference ratios (Eq.A.4.8-9).Bacteria take up LDOM in the way that the ratio of LDON (LDOP) to LDOC uptake is the same as the bulk N/C (P/C) ratio of the LDOM (Eq. A.4.15,Eq. A.4.21).Bacteria take up SDOM with higher N/P ratios to reflect that SDOM with higher N/P ratios is more labile (Eq. A.4.13).The ratio of SDON to SDOC uptake by bacteria would vary between the bulk N/C of SDOM and the bacterial reference cellular quota (Eq. A.4.16,A.4.22).Bacteria are modelled to either take up or release NH4 and PO4 to maintain their stable and consistent stoichiometry (Kirchman, 2000).Bacteria take up NO3 only if their cellular N/C ratio is smaller than their reference ratio (i.e., when bacteria are in short of nitrogen), in order to reflect higher energetic cost of NO3 uptake than NH4, but the amount of NO3 uptake is modelled to be no more than 10% of N-specific bulk L-and SDOM uptake, and the sum of NO3 and NH4 uptake is modelled to be no more than N-specific bulk L-and SDOM uptake (Eq.A.4.17-20).These limit the maximum NO3 uptake rate and set the inhibition of NH4 uptake on NO3 uptake.Bacteria excrete RDOM by transforming LDOM to .Bacteria also adjust their cellular stoichiometry by remineralizing NH4 and PO4 if carbon is in short (i.e., N and P in excess; and by excreting SDOC if carbon is in excess (i.e., nitrogen and phosphorus are in short; Eq.A.4.38-43).Bacteria are grazed by microzooplankton , and a certain percentage of bacteria gets lost to LDOC pool due to viral attack (Eq.A.4.47-49).

Physical forcing
The WAP-1D-VAR model v1.0 is forced by mixed layer depth (MLD), PAR at the ocean surface, sea-ice concentration, water-column temperature, vertical velocity, and vertical eddy diffusivity, at a temporal resolution of 1 day.Temperature, sea ice, and vertical eddy diffusivity are set up at every vertical grid (depth) point.
MLD is determined based on a finite difference density criterion with a threshold value of ∆σθ = 0.03 kg m -3 (Montégut et al., 2004) after calculating potential density of water mass from temperature and salinity Conductivity-Temperature-Density (CTD) observations.Vertical velocity, w, is assigned as zero because it is very weak in the surface waters of the study site and materials are transported vertically mostly by diffusion.The vertical eddy diffusivity scheme treats the rapid vertical mixing in the surface boundary layer by homogenizing model state variables instantaneously in the mixed layer (i.e., by averaging at every time step).Thus, Kz value above MLD is not required, and only Kz below MLD is calculated as follows: where z is depth (m) below MLD and KZ0 is the vertical eddy diffusivity at the bottom of the mixed layer (1.1 ´ 10 -4 m 2 s -1 ) (Klinck, 1998;Smith et al., 1999), and α is 0.01 (m -1 ).
Daily surface downward solar radiation flux (National Centers for Environmental Prediction reanalysis daily averages) is used to calculate sea surface PAR.PAR is estimated as 46% of the total solar radiation (Pinker & Laszlo 1992, Kirk 1994).The attenuation of PAR as a function of depth is calculated as follows: where z is depth (m), PAR0 is PAR level at sea surface (W m -2 ), kw is the attenuation coefficient for seawater (m -1 ), kc is the attenuation coefficient for Chl ((mg Chl) -1 m 2 ), and CHL is the Chl concentration (mg Chl m -3 ).
Sea-ice conditions in the coastal WAP do not necessarily represent solely local temperature and climate conditions, given that sea ice can be impacted by temperature, mixed layer, heat fluxes, regional winds, and other physical processes (Saenz et al. in review).We implement a sea-ice model scheme to account for light transmission through sea ice (5% of incident irradiance, as a typical transmittance value used in the Community Earth System Model) and non-linearities in the photosynthesis-irradiance (P-I) response under partial ice concentration (Long et al., 2015) using percent daily sea-ice concentration data (GSFC Bootstrap versions 2/3, derived from SMMR/SSMI satellite temperature brightness data binned into 25 by 25 km grid cells).In many previous models, the light-limitation term ℒ(I) is calculated as a function of mean irradiance I ̅ averaged over both ice-covered and open-water conditions, so ℒ(I ̅ ); instead we compute the mean of light-limitation term (ℒ(I) ''''' ) as a function of fractional sea ice and open water and incident irradiance: (5) where P C is the C-specific photosynthetic rate (d -1 ), P C MAX is the maximum photosynthetic rate (d -1 ), Ik is the parameter describing the light-saturation behaviour of the PI-curve (W m -2 ), Io is the open-water irradiance, Ii is the under-ice irradiance (i.e., Ii = 0.05 ´ Io), fi is the fraction of area covered with sea ice, and fo is the fraction of open water (i.e., fo = 1 -fi).

Variational data assimilation
The WAP-1D-VAR model v1.0 is equipped with a built-in data assimilation scheme based on a variational adjoint method (Lawson et al., 1995).This method generates optimal model solutions that minimize the difference between model results and observations by objectively optimizing model parameter values (Friedrichs, 2001;Spitz et al., 2001;Ward et al., 2010).In detail, the assimilation scheme (Figure 2) consists of four steps (Glover et al., 2011): 1) starting with initial values of the model parameters (see below), the model is integrated forward in time from specified initial conditions to calculate the difference between the model simulation and the field data, or the model-observation misfit (i.e., cost function; section 2.5, Eq. 6); 2) an adjoint model constructed using the Tangent linear and Adjoint Model Compiler (TAPENADE) is integrated backward in time to compute the gradients of the total cost with respect to the model parameters; 3) the computed gradients are then passed to a limited-memory quasi-Newton optimization software M1QN3 3.1 (Gilbert & Lemaréchal, 1989) to determine the direction and optimal step size by which the selected model parameters (see below) need to be modified in order to reduce the total cost; and 4) a new forward in time simulation is conducted using the new set of modified (optimized) parameter values.These four-step procedures are conducted in an iterative manner until the pre-set convergence criteria (i.e., low gradients of the total cost function with respect to optimized parameters and positive eigenvalues of the Hessian matrix) are satisfied to ensure that the optimized parameters converge and the total cost function reaches a local minimum.
Initial values of the model parameters (total of 72 free or optimizable parameters, Table 1) are assigned based on literature values (Caron et al., 2000, Luo et al., 2010, Garzio et al., 2013) without examining the effects of the initial parameter values on the model results prior to optimization.As is typical for many types of ecosystem models, a collection of what appear to be reasonable initial parameter estimates can result in relatively poor overall system behaviour because of system-level interactions of different model components.In most marine ecosystem models, these initial parameter values are subjectively and manually adjusted to improve the simulation, and the simulations with the initial, unadjusted parameter values are rarely shown.However, here with a more objective optimization approach that we conduct, the initial and optimized solutions can be explicitly compared (section 4).Optimization starts by submitting a subset of the 72 free model parameters rather than submitting all of them at once.This initial parameter subset consists of 10 different model parameters, the change of which yields the largest decrease in the total cost function (i.e., which also happens to be usually one per each state variable).These include αDA (initial slope of photosynthesis vs. irradiance curve of diatoms, mol C (g Chl a) -1 d -1 (W m -2 ) -1 ), αCR (initial slope of photosynthesis vs. irradiance curve of cryptophytes, mol C (g Chl a) -1 d -1 (W m -2 ) -1 ), Θ (maximum Chl/N ratio, g Chl a (mol N) -1 ), μBAC (maximum bacterial growth rate, d -1 ), r A max,BAC (maximum bacterial active respiration rate, d -1 ), gBAC (half-saturation density of bacteria in microzooplankton grazing, mmol C m -3 ), μMZ (maximum microzooplankton growth rate, d -1 ), μKR (maximum krill growth rate, d -1 ), and remvKR (krill removal rate by higher-trophic levels, (mmol C m -3 ) -1 d -1 ; Table 1).
When computed at the minimum of the cost function value, the inverse of the Hessian matrix provides the uncertainties of optimized parameters, cross-correlations among parameters, and sensitivities of the total cost function to each parameter (Matear, 1996;Tziperman & Thacker, 1989).High off-diagonal values in the inversed Hessian matrix indicate highly cross-correlated model parameters, so one of the highly cross-correlated parameters is removed from the optimization.The square root of a diagonal element in the inversed Hessian matrix is the logarithm of the relative uncertainty (σf) of the corresponding optimized parameter.The absolute uncertainty of the constrained parameter is calculated as pf ´ e ±σf where pf is the value of the optimized parameter (Table 1).If optimized to ecologically unrealistic values, parameters are kept back to their respective initial values and removed from the next optimization cycle.Optimized parameters with σf larger than 50% are updated but removed from the next optimization cycle (i.e., defined as 'optimized' parameters), while optimized parameters with σf smaller than 50% are updated and kept for the next optimization cycle (i.e., defined as 'constrained parameters').Constrained parameters are reported with the uncertainties, while optimized parameters are reported without the uncertainties (Table 1) because both changed parameters consist of an optimized model parameter set, but the parameters reported with the uncertainty ranges are the ones optimized with relatively small uncertainties and considered constrained.This way, a part of the initial parameter subset forms a final optimized parameter set.The gradients of the total cost function with respect to all 72 parameters are then evaluated, the parameters with large gradients (e.g., > 5) are re-submitted to optimization to further reduce the total cost, the gradients are evaluated again, and these cycles repeat until the termination of optimization.Optimization terminates when the gradients are reasonably low (e.g., < 10 -2 for constrained parameters, < 5 for optimized parameters, and < 10 for unoptimized parameters).This final optimized model parameter set forms the basis of the results presented throughout this study (section 4).Additionally, in order to assess the sensitivity of the model optimization results with regard to the initial parameter choice, we perturb by ± 50% a subset of the initial parameter values used in the reference (original) optimization experiments to form different initial parameter sets (a total of 15 sets consisting of partially or fully perturbed 18 parameters, Tables B1-2) and conduct new optimization experiments from each set (section 4.1).

Cost function
To represent a misfit between observations and model output, a total cost function is calculated as follows (Luo et al., 2010): where m and n represent assimilated data types and data points, respectively, M and Nm are the total number of assimilated data types and data points for data type m, respectively, σ m is the target error for data type m, am,n is observations, and a ( m,n is model output.Given the high biological productivity of the WAP waters and the approximate log-normal distribution of many marine biological variables, the base-10 logarithms of Chl and primary production (PP) are used in the cost function calculation to capture phytoplankton dynamics (Campbell, 1995;Glover et al., 2018).The target error is calculated for each data type as follows: ''''' is the climatological mean (over the select 9 growth-seasons, see below) of the observations and CV m is the averaged coefficient of variation (CV) of the observations of each data type in the mixed layer (due to observational error and seasonal and interannual variations) calculated using all of the observational data over 9 growth-season periods between 2002-2003 and 2011-2012, except the 2007-2008 growth season due to its missing data.These 9 growth seasons are chosen, instead of the multi-decadal observations available from Palmer LTER (since 1991), due to the relatively more complete data coverage in those seasons.The standard deviations are used as target errors of the log-converted data types.The CV of the log-converted data type is estimated as the average of ± 1 standard deviation in log space converted back into normal space (Doney et al., 2003;Glover et al., 2018).Hereafter, we present the total cost normalized by M (J equivalent to J/M hereafter) as it indicates the model-observation misfit equivalent to a reduced Chi-square estimate of model goodness of fit.We report the normalized total cost J along with normalized costs of individual data types throughout this study.J = 1 indicates a good fit, J >>1 indicates a poor fit or underestimation of the error variance, and J <<1 indicates an overfitting of the data, fitting the noise, or overestimation of the error variance.

Modelling framework
To examine the applicability of the WAP-1D-VAR model v1.0 to the coastal WAP region, we select a nearshore Palmer LTER water-column time-series station, Station E (64.77°S, 64.05°W), as the modelling site that is ~200 m deep and situated approximately 3 km south of Palmer Station and 6.5 km northeast of the head of Palmer Deep (Sherrell et al., 2018).Physical forcing (Figure 3) and data types assimilated are derived from roughly semi-weekly physical, chemical and biological profiles collected from small boat via a profiling CTD and discrete water samples at Station E. When weather and ice conditions permit, water column sampling at the station has been conducted twice a week over the growth season.Seven upper-ocean layer depths (2.5, 10, 20, 30, 40, 50 and 60 m) are chosen for the model vertical grids.The model depth can be extended to as deep as needed, but this study is focused to upper 60-m water column to fully take advantage of the large data availability.Also, conceptually, the application of the 1-D model framework makes the most sense for the upper water column dominated by local seasonal processes, and extension of the model into deeper water well below the maximum seasonal mixed layer becomes more problematic because of the growing importance of lateral advective processes that are not well captured in the 1-D model framework.The vertical structure of the water column can be affected by growing sea ice due to reduced winddriven turbulence and brine rejection during winter, but this is what a prognostic, coupled ocean-ice 1-D model can offer to simulate, not our diagnostic forcing based model used in this study.Also, because our model simulates only the spring-summer growth season, the impact of winter sea ice on ecosystem dynamics is less of a concern.
Given the routine observations of Palmer LTER available over the growth season (October -March), we simulate one example growth season with the most complete data coverage, from October 2002to March 2003(2002-2003 growth season hereafter), instead of a series of different growth seasons in a continuous manner.The example growth season simulations utilize this year's specific observed physical forcing fields and assimilated biological and biogeochemical observations.Each Palmer LTER growth season should be modelled to have its own unique optimized parameter set, as well as initial conditions and physical forcing that together determine the model solution for that year; however, only 2002-2003 growth season simulations are modelled in this study for model analysis and evaluation.

Initial and boundary conditions
Model initial conditions are prescribed 135 days before the model start date for the growth season (October 15, 2002), so on June 1, 2002.This 135-day spin up is conducted to minimize the impact of initial conditions on the model output over the growth season.Initial conditions are prepared by optimizing the full growth seasonal cycle forced by climatological physics and assimilated with climatological observations and with the same bottom boundary conditions used in the optimization of the 2002-2003 growth season (i.e., climatological model; using climatological physics and observations averaged over 9 growth-season periods between 2002-2003 and 2011-2012 except the 2007-2008 growth season due to its missing data).For the first climatological model simulation initial conditions are prepared by adjusting manually following literature values (e.g., Luo et al. 2010).Due to strong interannual variability in the phytoplankton bloom phenology at Palmer Station, averaging across all these 9 years does not reflect distinct seasonal phytoplankton peaks, leading to underestimated phytoplankton values (not shown).To capture this non-linear aspect of the coastal WAP system, we construct the climatological year by applying a single time shift to all variables so that a seasonal PP peak of each year lines up with an average date of seasonal PP peaks from all years.Most biological initial conditions on June 1 are close to zero given the lack of active physiological processes in the very low light and the presence of sea ice during wintertime before the model growth season starts.All the data types are set to zero at the lower boundary (bottom) except for NO3, PO4, SDOC, SDON, and SDOP in which the climatological values at 65 m are used for lower boundary values (25.9 mmol m -3 , 1.9 mmol m -3 , 6.5 mmol m -3 , 0.6 mmol m -3 , and 0.03 mmol m -3 , respectively).

Assimilated data
We include the data types directly related to corresponding model outputs, including a mix of ecosystem stocks or state variables -NO3, PO4, Chl for diatoms and cryptophytes, bacterial biomass, microzooplankton biomass, SDOC, POC, and PON as well as carbon flows among model stocks -bulk net PP and bacterial production (BP).These data sets have been sampled semi-weekly at Palmer Station E (64.77°S, 64.05°W), the same location where our model is set up, and are available from the Palmer LTER data website (see Data availability).The distinction between diatoms and cryptophytes is established by assimilating phytoplankton taxonomic-specific Chl data for diatoms and non-diatom species derived from a High-Performance Liquid Chromatography (HPLC) and CHEMTAX analysis (Schofield et al., 2017), but given cryptophytes being the second dominant species in the water samples at the study site, cryptophytes are assumed to represent all non-diatom species for modelling purposes.Given that POC (PON) from bottle filtration may capture both living biomass and detrital material, we adjust the observed POC (PON) by subtracting phytoplankton and bacterial C (N) biomass to estimate the detrital pool, in order to only include non-living particles to detrital pool.When phytoplankton or bacterial biomass data are not available, we assign climatological (2002-2003 to 2011-2012) fractions of POC (PON) to detrital pool.Phytoplankton-and bacterial biomass accounts for 26% of total POC and 29% of total PON.In converting Chl to phytoplankton carbon (nitrogen) biomass, the maximum Chl/C (Chl/N) ratio submitted for optimization is used along with other reference ratios (Table 1).Microzooplankton biomass data are not available for the full time-series, so their data from grazing experiments at Palmer Station (Garzio et al., 2013) are assimilated to at least provide constraints on bacterial and cryptophyte grazing processes.However, due to the discrepancy in the timing and location from model simulations of this study, the microzooplankton modelobservation misfits are not analysed in the present study.Krill biomass data are not assimilated due to the strong patchiness of their distribution that may hinder proper model optimization.The vertical profiles of most of the data types are assimilated, whereas average NO3 and PO4 concentrations in the mixed layer are assimilated due to the difficulty of simulating depthdependent nutrient concentrations and the fact that net PP is mostly determined by surface nutrient concentrations (Luo et al., 2010).BP (mmol C m -3 d -1 ) is derived from the 3 H-leucine incorporation rate (pmol l -1 h -1 ) data using the conversion factor of 1.5 kg C (mol leucine) -1 incorporated (Ducklow, 2000).Bacterial biomass (mmol C m -3 ) is estimated from bacterial abundance measured by flow cytometry with the conversion factor of 10 fg C cell -1 (Fukuda et al., 1998).SDOC is calculated by subtracting the background concentration (41.2 mmol m -3 for the modelling site) from total DOC concentration.

Uncertainty analysis
Uncertainties of the optimized parameters are computed from a finite difference approximation of the complete Hessian matrix (i.e., second derivatives of the cost function with respect to the model parameters) during the iterative optimization process (details in section 2.4).We then conduct Monte Carlo experiments to calculate the impact of the optimized parameter uncertainties on the model results.We first create an ensemble of parameter sets (n = 1000) by randomly sampling values within the uncertainty ranges of the constrained parameters, and then perform a model simulation using each parameter set.1000 Monte Carlo experiments were shown to be adequate from a series of tests with different numbers of Monte Carlo sampling (n = 500-2000), where standard deviations of model simulated values converged after >1000 Monte Carlo sampling (not shown).All uncertainty estimates are calculated following standard error propagation rules and presented as ± 1 standard deviation in the study.

Model skill assessment
In the case of the example growth season (2002)(2003) modelled in this study, the iterative data assimilation-parameter optimization procedure reduced by 58% the misfits between observations and model output compared to the misfits obtained from the initial parameter values (Table 2).The optimized model solution satisfied the pre-set convergence criteria, with the low gradients of the total cost with respect to the optimized parameters and positive eigenvalues of the Hessian matrix.Notably, this was achieved by optimizing a subset of 12 (9 constrained and 3 optimized) parameters among the total of 72 optimizable parameters (Table 1, section 4.2).To examine the sensitivity of the optimized model solution to the initial parameter choice, a series of new optimization experiments (n = 15) were conducted with a varying subset of the initial parameter values perturbed by ± 50% of those used in the original optimization experiment (Table B1).These experiments showed that the optimized model results (i.e., the reference case; Table 1) were not sensitive to the initial choice of the parameters.The 15 different initial parameter sets resulted in a range of initial model-observation misfits, some substantially larger than the reference case (14.25-28.24versus 14.85 for the reference case).However, the total normalized optimized cost values of the 15 sensitivity experiments (5.79-7.19)were similar to that of the reference case of 6.42.In the sensitivity experiment #12, the initial modelobservation misfit was ~2 times larger than that of the reference case, and there was up to 76% of the reduction in the modelobservational misfit (versus 58% of the reduction in the reference case; Table B1, Table 1).These results suggest that no matter where in parameter space the optimization process starts from, the optimization scheme takes the model cost function to similar local minima.Importantly, this was achieved by similar subsets of the optimized parameters: μDA, μCR, r A max,BAC, and μMZ were optimized in all cases, while αDA, αCR, Θ, μBAC, gCR, g'DA, gBAC, μKR, gMZ, and remvKR were optimized except for a few cases (Table B1).The uncertainties of the optimized parameters were also similar among different optimizations, with most of the relative errors < 0.5.Constrained parameter values and their uncertainty ranges averaged over the sensitivity experiments (Table B2) were comparable to those in the original optimization experiments (Table 1).
Overall, there was a good model-data fit with the largely decreased cost value for each data type after optimization (Table 2).Optimization yielded Jf close to 1 for all data types, compared to the initial model solution where three data types -diatom Chl, crypto Chl, and bacterial biomass -had particularly poor model fits to observations and underestimated error variances (J >> 1).Compared to the initial (unoptimized) model results, the average errors (εbias, Doney et al., 2009;Stow et al., 2009) in the optimized model results indicated that diatom Chl, cryptophyte Chl, bacterial biomass, BP, and POC had reduced model biases, while NO3, PO4, PP, SDOC, and PON had increased model biases (for both positive and negative biases, defined as εbias > 0 and εbias < 0, for model overestimation and underestimation of the observation, respectively).Optimization resulted in the negative model bias for PO4, compared to the positive model bias in the initial model results.The point-to-point comparison plots showed that there were seasonally consistent, negative model biases for PP, POC, and PON (Figure B1).Model skill was further evaluated with a Taylor diagram (Taylor, 2001) summarizing the statistics of the correlation coefficient between model output and observations, normalized standard deviation (by the standard deviation of the observations), and centred (bias removed) root-mean-square difference (RMSD) for each data type, in which a better model skill is characterized by a higher correlation, a normalized standard deviation close to 1, and a lower RMSD (Figure 4).Optimization resulted in better model skills for cryptophyte Chl, PP, BP, and bacterial biomass via increased correlation coefficients and lowered RMSD (Figure 4B), compared to those in the unoptimized model results (Figure 4A).After optimization the normalized standard deviations of PP, BP, bacterial biomass, phosphate, POC, and PON were closer to 1 (Figure 4B).Direct comparisons with the observational data showed that the optimized model parameter set captured better the increases in diatom biomass early in the season, cryptophyte biomass in January, and bacterial biomass in mid-February, compared to the unoptimized model parameter set (Figures 5A-B, Figure B2).

Optimized parameters
The number of the optimized parameters in this study is small and comparable to those from other data-assimilative model focused on different marine environments (Friedrichs, 2001;Friedrichs et al., 2006Friedrichs et al., , 2007;;Luo et al., 2010).This is consistent with the general behaviour of marine plankton ecosystem models, in which well-posed model solutions would be found with only a subset of independent model parameters due to many cross-correlated parameters inherent in nonlinear model equations (Fennel et al., 2001;Harmon & Challenor, 1997;Matear, 1996;Prunet et al., 1996).Ecosystem models with a relatively large number of unconstrained parameters (i.e., equivalent to the optimized parameters with high uncertainties in the present study) might reduce total costs to a greater extent, but could possess low predictive skill as a result of being overtuned to fit noise in the observations (Friedrichs et al., 2007).Also, there are several field and lab-based studies at the study site or in a similar polar environment that reported the values of the model parameters used in the WAP-1D-VAR model v1.0, including the bacterial growth rate of 0.82 d -1 , total phytoplankton (including large cells like diatoms) growth rate of 0.33-0.55d -1 , nanophytoplankton (corresponding to cryptophytes) growth rate of 0.52-0.99d -1 (Garzio et al., 2013), and the microzooplankton growth rate of up to 1.0 d -1 (Caron et al., 2000).The optimized values of the maximum bacterial, diatom, cryptophyte, and microzooplankton growth rates in this study were 1.06 d -1 (0.93-1.20 d -1 ), 0.77 d -1 (0.68-0.88 d -1 ), 0.72 d -1 (0.61-0.85 d -1 ), and 1.18 d -1 (1.10-1.26d -1 ), falling in the ranges of those measured bacterial, total phytoplankton, nanophytoplankton, and microzooplankton growth rates, respectively.
The ensemble of the model state variables and flows obtained from the Monte Carlo experiments had generally small standard variations at each model time step and grid, suggesting the robustness of the modelled fields against the variations in the optimized parameter values (Figures B3, B4).To examine the translation of the optimized parameter values to altered functioning of the WAP biogeochemical processes, we compared two different sets of the model simulation results -one based on the initial parameter values (Figures 5A, 6A, 7A) and the other based on the optimized parameter values (Figures 5B, 6B,  7B).However, due to the non-linearities in the model it is not straightforward to identify what causes the parameter variations, except for a few cases in which the changes in the parameter values are clearly linked to the difference in the model state variables and flows.The first case is the relation of the increased gBAC value (bacterial half-saturation concentration in microzooplankton grazing, mmol C m -3 ) to the elevated bacterial accumulations after optimization (Table 1, Figures 5, 7).The second case is the link between Θ (maximum Chl/N ratio, g Chl a (mol N) -1 ) and the relative dominance of cryptophytes in total phytoplankton accumulations.It has been demonstrated that the variations of Θ are driven by an imbalance between the rate of light absorption and energy demands for photosynthesis and biosynthesis in phytoplankton cells (Geider et al., 1997).Θ can also change because of the variations in phytoplankton photo-acclimation or physiological differences across phytoplankton groups, from a lower Θ value for smaller species to a higher Θ value for larger diatom cells (Geider, 1987).Θ was optimized to a 30% lower value than the initial parameter value (Table 1), in order to simulate the relatively larger proportion of cryptophytes in total phytoplankton accumulations in the optimized model results compared to the unoptimized model results (Figures 5, 7).By contrast, the remaining cases are not as clear because the first-order impact of parameter variations on the model results is less direct and more nuanced.Compared to the unoptimized results, the decreases in μDA (diatom C-specific maximum growth rate, d -1 ), μCR (cryptophytes C-specific maximum growth rate, d -1 ), αDA (initial slope of P-I curve of diatoms, mol C (g Chl) -1 d -1 (W m -2 ) -1 ), and αCR (initial slope of P-I curve of cryptophytes, mol C (g Chl) -1 d -1 (W m -2 ) -1 ) did not lead to decreased diatom and cryptophyte accumulations, presumably due to decreased gMZ (microzooplankton half-saturation concentration in krill grazing, mmol C m -3 ) and increased remvKR (krill removal rate by higher-trophic levels, (mmol C m -3 ) -1 d -1 ) after optimization (Table 1, Figures 5, 7).Similarly, the decreased μBAC (maximum bacterial growth rate, d -1 ) and the increased r A MAX,BAC (bacterial maximum active respiration rate, d -1 ) did not lead to decreased bacterial accumulations, presumably due to the increased gBAC (bacterial half-saturation concentration in microzooplankton grazing, mmol C m -3 ) and the decreased gMZ (microzooplankton half-saturation concentration in krill grazing, mmol C m -3 ).

Ecosystem indices
We calculated key ecosystem indices for the modelled growth season, including NPP (directly comparable to 14 C-PP observations), net community production (NCP; i.e., NCP = NPP -bacterial-, microzooplankton-, and krill respiration), BP, and POC export (sinking) flux (Figure 6).Setting an upper limit for lateral or vertical carbon export from the euphotic zone (Dugdale & Goering, 1967), over appropriate time and space scales NCP is quantitatively equivalent to new production that is supported via external sources of nitrogen (Ducklow & Doney, 2013).In both optimized and unoptimized model results, NPP increased after complete sea-ice retreat, but a brief ice-edge bloom was simulated under sea ice at the beginning of the growth season (Figures 3, 6).Seasonal patterns of NCP resembled those of NPP and occasionally fell below zero (i.e., the net heterotrophy) in subsurface waters for both optimized and unoptimized cases (Figure 6).The POC export flux increased over time and reached the maximum value at the end of the growth season in both model results, but there were two major POC flux events separated by weaker, in-between flux events in December in the optimized results that the initial model results did not capture (Figure 6b).After optimization, the correlation coefficients adjusted from 0.88 to 0.36, 0.89 to 0.68, and 0.45 to 0.73 for the NPP-vs.-NCPpair, the NPP-vs.-BPpair, and the NCP-vs.-POCexport flux (lagged by 30 days) pair (all p < 0.001).In the optimized model results, the growth-season mean of the depth-integrated NPP, NCP, and BP in the 60 m water column, and the 30-day lagged POC export flux at 60 m were 19 , and 2 ± 0.3 mmol C m -2 d -1 (uncertainties propagated from season-averaging in Figure 6b and Monte Carlo uncertainties in Figure B4), compared to 28 , and 2 ± 0.2 mmol C m -2 d -1 in the unoptimized model results (uncertainties from season-averaging in Figure 6a).
The mean e-ratio, defined as the growth-season mean of the 30-day lagged POC export flux divided by the growthseason mean NPP (i.e., particle export efficiency), was 0.11 ± 0.05 (uncertainties propagated from season-averaging in Figure 6b and Monte Carlo uncertainties in Figure B4) in the optimized model results, compared to 0.07 ± 0.02 (uncertainties from season-averaging in Figure 6a) in the unoptimized model results.The mean f-ratio, defined as the amount of NO3 uptake divided by the amount of NO3 and NH4 uptake both, was 0.88 ± 1.52 in the optimized model results, compared to 0.84 ± 0.19 in the unoptimized model results (not shown).The higher mean f-ratio relative to the mean e-ratio in this study implies an imbalance between production and export at the study site, at least during the modelled period.Excess new production relative to export production (as derived from sediment traps and 234 Th disequilibrium; Ducklow et al., 2018) was previously observed in the WAP, presumably due to diel vertical migration, DOM export, lateral export, and diffusive loss of PON via diapycnal mixing (Stukel et al., 2015).Stukel et al (2015) reported up to 5 times larger new production via NO3 uptake than export production via Th-based N export along the coastal WAP.Several additional mechanisms might be responsible for driving the discrepancy between production and export.First, given that the assimilated pool of suspended POC in the model formulation is not a good indicator of a rapidly sinking detrital pool dominating particle export, the WAP-1D-VAR model v1.0 does not capture large, short-lived particle flux events (e.g., fecal pellets produced by a large swarm of krill), underestimating POC export flux.Second, the WAP-1D-VAR model v1.0 export scheme does not consider DOC export that would lower the production-export discrepancy.Finally, RDOC is not explicitly modelled in the model, due to its much longer time scale than the model time scale, so accumulated and not-exportable RDOC pool would contribute to the deviation of the modelled e-ratio from the modelled f-ratio.Indeed, the modelled mean e-ratios in our study, for both optimized and unoptimized cases, situate at the lower end of the range of the e-ratios measured or estimated in the WAP waters (Ducklow et al., 2018;Sailley et al., 2013;Stukel et al., 2015;Weston et al., 2013), but optimization increased the e-ratio by 60% and thus made it closer to the literature values.
The mean BP/NPP ratio was 0.05 ± 0.06 (uncertainties propagated from season-averaging in Figure 6b and Monte Carlo uncertainties in Figure B4) in the optimized model results, compared to 0.11 ± 0.04 (uncertainties from season-averaging in Figure 6a) in the unoptimized model results.The modelled mean BP/NPP ratio for both optimized and unoptimized cases correspond well to the estimates from other measurement-and observation-based studies (Ducklow et al., 2012;Kim & Ducklow, 2016).Relatively low bacterial activity in productive Antarctic waters, typically reflected as a low BP/PP ratio, has been attributed to low LDOM availability for bacterial growth (Kirchman et al., 2009), low temperature (Pomeroy & Wiebe, 2001), or top-down control via grazing and viral lysis (Bird & Karl, 1999).

Mean carbon stocks and flows
We summarized the growth-season means of the carbon stocks and flows in the entire food web (Figure 7).The WAP-1D-VAR model v1.0 captured several key WAP ecological and biogeochemical features, including the lack of macronutrient limitation (NO3 and PO4 drawdown by phytoplankton utilization but remaining well above their half-saturation constants, Table 1) and comparable values of the assimilated and non-assimilated model state variables (Ducklow et al., 2007(Ducklow et al., , 2012(Ducklow et al., , 2018;;Kim et al., 2016;Moline et al., 2008;Smith et al., 2008), providing confidence in the model simulations.For instance, growth-season measurements in 2017-2018 at Palmer Station showed a strongly patchy krill distribution, with the mean biomass of 0.12 ± 0.04 mmol C m -3 and the maximum biomass of 0.57 mmol C m -3 when krill were present (unpublished data provided by D. Steinberg), falling in the range of the modelled krill biomass values (0.13 ± 0.03 mmol C m -3 ; calculated from Figure 7b).The WAP-1D-VAR model v1.0 also simulated several important ecosystem metrics comparable to other statistical modelling studies.For instance, the modelled phytoplankton seasonal patterns in the present study are consistent with physicochemical attributes revealed by a distinct ecological seascape pattern in the coastal WAP (Bowman et al., 2018), including low Chl and high nutrients in the first half of the growth season followed by high Chl and low nutrients in the second half of the growth season.A steady-state solution based, inverse modelling study quantified different food-web states using ecosystem network indices from Palmer LTER annual summer cruises along the WAP shelf region (Sailley et al., 2013).Their network indices include the ratio of C export to total PP (i.e., equivalent to e-ratio in our study) and the ratio of recycling (the sum of flows into respiration and DOC pool) to total PP, where more (less) recycling favourable microbial food-webs are characterized by greater (smaller) ratios of recycling to total PP and smaller (greater) ratios of total C export to total PP (Legendre & Rassoulzadegan, 1996).As discussed above, the modelled mean e-ratio in the present study is smaller than the estimates in the inverse modelling study for the WAP shelf region (Sailley et al., 2013), but consistent with their conclusion on the food-web status of the modelled growth season (2002)(2003) positioned on the microbial food-web side.The discrepancy in the e-ratio values between the present study and Sailley et al. (2013) may be attributed to fundamentally different model formulation (i.e., time-evolving modelling for the WAP-1D-VAR model v1.0 versus steady-state modelling) and optimization approach, or due to relatively strong microbial food-web activity at our coastal site compared to the shelf region.Microbial food-web activity can be approximated by quantifying the amount of fixed carbon flowing through heterotrophic bacteria (Carlson et al., 1999;del Giorgio & Cole, 1998;Ducklow, 2000;Ducklow et al., 2012).According to this approach, microbial food-web activity from the optimized model results was around 38 ± 16%, calculated as the ratio of bacterial L-and SDOC uptake to PP (i.e., (arrow 13 + arrow 14)/arrow 1 in Figure 1, mean ± uncertainties from season-averaging and Monte Carlo uncertainties in Figure 7B).On average, SDOC supported 1 ± 2% of the total bacterial C uptake, or C demand (i.e., arrow 14/(arrow 13 + arrow 14) in Figure 1, mean ± uncertainties from season-averaging and Monte Carlo uncertainties in Figure 7b), but could be an important bacterial C source when LDOC became scarce as the growth season progressed (Figure 5b).Indeed, several observational studies speculated that the WAP bacteria utilize SDOM in short of LDOM (Ducklow et al., 2011;Kim & Ducklow, 2016;Luria et al., 2017).

Summary
We developed the WAP-1D-VAR model v1.0, a one-dimensional variational data assimilation model specific to the coastal WAP region, evaluated the model performance and robustness using a variety of quantitative metrics, and discussed the model applicability with regard to capturing the key WAP ecological and biogeochemical features using the data from the example growth season.The data assimilation scheme significantly reduced the model-observation misfits via the optimized model parameter set that adjusted the simulation to better match the Palmer LTER observations in 2002-2003.We also explored the nuanced question of how the observations influenced the data assimilation process, drove the variations in optimized parameter values relative to their corresponding initial parameter values, and affected the resulting model simulations.The WAP-1D-VAR model v1.0 successfully simulated the variables and flows not included in data assimilation, with the values consistent and comparable with other measurement-and observation-based studies in the coastal WAP.Importantly, the data assimilation scheme enabled the available observational data to constrain processes that were poorly understood, including the partitioning of NPP by different phytoplankton groups, the optimal Chl/C ratio of the WAP phytoplankton community, and the partitioning of DOC pools with different lability.Up to this point, a range of observational studies has provided snapshots of ecosystem and biogeochemical processes in the WAP.Yet, we have little understanding of the driving processes that underlie the connections between each component in complex food-web interactions.We used dataassimilative modelling to glue these snapshots together to explain better the observed dynamics and further understand the previously poorly constrained processes in the coastal WAP system., d -1 1.00 0.72 (0.61-0.85) 2.51´10 -4 αDA, Initial slope of P-I curve of diatoms, mol C (g Chl) -1 d -1 (W m -2 ) -1 0.30 0.13 (0.10-0.19) -1.55´10 -4 αCR, Initial slope of P-I curve of crypto., mol C (g Chl) -1 d -1 (W m -2 ) -1 0.20 3.89´10 -2 (n/a) 0.45 βDA, Light inhibition parameter for diatom photosynthesis (W m -2 ) -1 5.00´10 -3 --1.10 βCR, Light inhibition parameter for crypto.photosynthesis (W m -2 ) -1 5.00´10 9.00´10 -4 --5.17´10 -2 q C N,RDOM, RDOM N/C ratio, mol N (mol C) -1 0.05 --0.10 q C P,RDOM, RDOM P/C ratio, mol P (mol C) -1 6.50´10 -4 -4.70´10 -3 q C N,POM, N/C ratio for POM production by microzoo.and krill, mol N (mol C) -1 0.12 -0.12 q C P,POM, P/C ratio for POM production by microzoo.and krill, mol P (mol C) -1 4.50´10 -3 -6.24´10 -2 rntrf, Nitrification rate (NH4 to NO3), d -1 7.60´10 N,MIN,DA, Minimum N/C ratio of diatoms 3.40´10 -2 q C N,MAX,DA, Maximum N/C ratio of diatoms 0.17 q C N,RDF,DA, Reference (Redfield) N/C ratio of diatoms 0.15 q C P,MIN,DA, Minimum P/C ratio of diatoms 1.90´10 -3 q C P,MAX,DA, Maximum P/C ratio of diatoms 1.59´10 -2 q C P,RDF,DA, Reference (Redfield) P/C ratio of diatoms 9.40´10 -3 q C N,MIN,CR, Minimum N/C ratio of crypto.3.40´10 -2 q C N,MAX,CR, Maximum N/C ratio of crypto.0.17 q C N,RDF,CR, Reference (Redfield) N/C ratio of crypto.0.15 q C P,MIN,CR, Minimum P/C ratio of crypto.1.90´10 -3 q C P,MAX,CR, Maximum P/C ratio of crypto.1.59´10 -2 q C P,RDF,CR, Reference (Redfield) P/C ratio of crypto.9.40´10 -3 q C N,BAC, Reference (optimal) N/C ratio of bacteria 0.18 q C P,BAC, Reference (optimal) P/C ratio of bacteria 0.02 q C N,MZ, Reference (optimal) N/C ratio of microzoo.0.20 q C P,MZ, Reference (optimal) P/C ratio of microzoo.2.20´10 -2 q C N,KR, Reference (optimal) N/C ratio of krill 0.20 q C P,KR, Reference (optimal) P/C ratio of krill 8.00´10 -3 ϵDA, Grazing limit to the amount of diatoms available for microzoo.grazing, mmol C m -3 1.00´10 -3 ϵCR, Grazing limit to the amount of crypto.available for microzoo.grazing, mmol C m -3 2.95    μDA, μCR, αDA, αCR, Θ, μBAC, gBAC, μMZ, μKR, gMZ  First, for each sensitivity experiment lower and upper bounds of the constrained parameter are calculated as pf ´ e +σf and pf ´ e -σf (where pf is the value of the constrained parameter and σf is the square roots of diagonal elements of the inverse of the Hessian matrix), respectively.Then we form the "lower (upper) bound parameter set" that only consists of the lower (higher) bounds of the constrained parameters from each experiment, and average those across the sensitivity experiments (n = 15) to calculate the lower (upper) bound listed in parentheses.1.00 1.00 (0.27) 0.15 (n/a) 0.34 (0.28-0.43) remvKR, Krill removal rate by higher-trophic levels, (mmol C m -3 ) -1 d -1 0.10 0.10 (0.02) 0.43 (n/a) 0.14 (0.14-0.15) wnsv, Detritus vertical sinking velocity, m d -1 5.00 4.83 (0.65) -diss, Detrital dissolution rate, d -1 0.14 0.14 (0.02) --

Figure 1 .
Figure 1.Ecosystem model.The WAP-1D-VAR model v1.0 is forced by five different physical forcing, denoted as a horizontal row across the top of the schematic.The ecosystem component incorporates eleven different prognostic state variables.Higher level and refractory dissolved organic matter (RDOM) are represented implicitly.605

Figure 2 .Figure 3 .
Figure2.Variational data assimilation.A variational adjoint scheme is employed for the parameter optimization and data assimilation processes (adapted fromGlover et al., 2011).Gradient: the sensitivity of the total cost function with respect to model parameter from optimization.610

Figure 4 .Figure 5 .Figure 6 .
Figure 4. Model skill assessment.The Taylor diagrams using a polar-coordinate system summarizing the model-observational 620 correspondence for each model stock and flow for the modelled growth season 2002-2003 before (a) and after optimization (b).The angular coordinate for the Pearson correlation coefficient (r), the distance from the origin for the standard deviation normalized by the standard deviation of the observation, and the distance from point (1,0), marked as REF on x-axis, for the centred (bias removed) root-mean-square difference (RMSD) between model results and observations.625

Table 1 . Summary of model parameters.
Summary of the model parameter symbol and definition, initial parameter values (p0) and optimized values (pf) for optimizable parameters, the cost function gradient with regard to the optimized parameter (∂J/∂p), and prescribed values for fixed model parameters over the course of simulations.The parameter with 'n/a' in the parenthesis is an optimized parameter with a large relative uncertainty, while the parameter with values in the parenthesis is a 655 constrained parameter (optimized with a low relative uncertainty) with its upper and lower bounds.The uncertainties for these upper and lower bounds are calculated as: pf ´ e ±σf where pf is the value of the constrained parameter and σf is the square roots of diagonal elements of the inverse of the Hessian matrix.The cost function gradient with regard to the optimized parameter (∂J/∂p) after data assimilation is defined as: ΔJ/e Δp where e Δp ≈ Δp for an infinitely small Δp.For example, a 10% change of a parameter (Δp = 10%) leads to a total cost change equivalent to 10% of the corresponding gradient.660 Conversion rate of SDOM to RDOM, d -1

Table 2 . Data types, observed means, coefficient of variation, target errors, and costs before and after optimization.
Stow et al. (2009)( '), coefficient of variation (CV), and target error (s) of each assimilated data type used for calculating 665 the cost function before and after optimization.J0 is the normalized cost function before optimization and Jf is the normalized cost function after optimization (Eq.6).Data type units: mmol m -3 for NO3, PO4, diatom Chl, cryptophyte Chl, bacterial biomass, SDOC, and POC; mmol N m -3 for PON; and mmol C m -3 d -1 for PP and BP.The average error (εbias) of each data type (for non-transformed or raw Chl and PP) is calculated fromStow et al. (2009)before and after optimization where a positive value indicates the model overestimation of the observation and vice versa.670

Table B2 . Summary of model parameters from initial parameter perturbation experiments.
Summary of the model parametersymbol and definition, initial parameter values (p0, ref)and optimized values(pf, ref)for the original (reference) experiments (Table1) and initial parameter values (p0) and optimized values (pf) averaged from the sensitivity experiments (n = 15).p0, ref and the mean p0 are the same for most parameters because of the perturbations by ±50% their original initial parameter values (with standard deviation in parentheses), while p0, ref and the mean p0 are different for AE and wnsv because of the perturbations only by -50% their original initial parameter values (with standard deviation in parentheses).Numbers in parentheses for pf are the uncertainty ranges (lower and upper bounds) averaged across the sensitivity experiments as follows.