GROWTH OF SACCHARINA LATISSIMA IN AQUACULTURE SYSTEMS MODELED USING DYNAMIC ENERGY BUDGET THEORY

Aquaculture is an industry with the capacity for further growth that can sustainably feed an increasing human population. Sugar kelp (Saccharina latissima) is of particular interest for farmers as a fast-growing species that benefits ecosystems. However, as a new industry in the U.S., farmers interested in growing S. latissima lack data on growth dynamics. To address this gap, we calibrated a Dynamic Energy Budget (DEB) model to data from the literature and a 2-year growth experiment in Rhode Island (U.S.). Environmental variables forcing model dynamics included temperature, irradiance, dissolved inorganic carbon (DIC) concentration, and nitrate and ammonium concentration. The modeled final estimate for S. latissima blade length (cm) was reasonably accurate despite underestimation of early season growth. Carbon limited winter growth due to a low modeled specific relaxation rate (i.e. the lightdependent reactions of photosynthesis) for some model runs; other model runs displayed nitrogen limitation which occasionally led to length overestimation and underestimation due to the degree of interpolation necessary from the field data. The model usage, however, is restricted to S. latissima grown in an aquaculture setting because of assumptions made about tissue loss, summer growth patterns, and reproduction. The results indicate that our mechanistic model for S. latissima captures growth dynamics and blade length at the time of harvest, thus it could be used for spatial predictions of kelp aquaculture production across a range of environmental conditions. The model could be a particularly useful tool for further development of sustainable ocean food production systems in the U.S. involving seaweed.


Introduction
With a growing global population, one of the greatest challenges is providing healthy diets from sustainable food systems (Duarte et al., 2009, Merino et al., 2012, Willett et al. 2019. Meeting these demands requires strategies for optimization and development across each food production industry in both terrestrial and aquatic habitats. Aquaculture is currently the fastest growing food production sector in the world and now produces more seafood than wild-capture fisheries (FAO, 2018).
Within aquaculture, the ocean is increasingly viewed as an important area of expansion. This is because mariculture, or the farming of seafood in the ocean, has had one of the largest relative production increases in the last thirty years and has positive growth trajectories in many countries (Cottrell et al., 2018;Gentry et al., 2019). In fact, mariculture systems produced around 28.7 million tons of food in 2016 (FAO, 2018). The species, methods, and location of aquaculture will dictate its future contribution to the global food system within sustainable limits, particularly in the ocean (Gentry et al. 2017, Willett et al. 2019).
Aquaculture can have negative environmental impacts. In open systems of fed species, this negative impact is largely due to concentrated flows of feces and feed wastage leading to eutrophication (Wu, 1995) and alteration of food webs (Herbeck et al., 2013). Food for aquaculture species can put more pressure on terrestrial or aquatic food production systems, where the biomass harvested of one species is then reconstituted and serves as the input for the farming of another species (Cottrell et al., 2018). Open aquaculture systems composed of species that do not require supplemental feed or nutrients (i.e., primary producers and filter feeders) avoid these harms and instead can perform important ecosystem services such as removing dissolved organic and inorganic nutrients (Alleway et al., 2019). In particular, seaweed is of interest for nutrient capture in integrated multi-trophic aquaculture (IMTA) systems where they utilize the inorganic carbon and nitrogenous compound outputs from fed species such as finfish (Chopin, 2015). Other potential benefits of growing seaweed are its ability to combat hypoxia from terrestrial food production systems and even protect shorelines through dampening of wave energy (Duarte et al., 2017).
Outside of these ecosystem services, growing seaweed has been proposed as a way to engage a wider public audience with climate change via offsetting carbon emissions (Froehlich et al., 2019). Seaweed aquaculture has the potential to generate net positive environmental and social impacts, but this industry has been traditionally concentrated in Asian countries (FAO, 2018).
Seaweed aquaculture in the U.S. is a nascent industry. The U.S. does not produce enough aquatic plants to even register in the global production statistics (< 0.1%; FAO, 2018). In the Northeast U.S., sugar kelp (Saccharina latissima) is a local species of recent interest for food, biofuel, bioremediation, and pharmaceutical products (Forbord et al., 2012). S. latissima grows to large sizes very quickly and individuals can reach lengths greater than four meters in less than five years in the wild (Borum et al., 2002). In a single season of aquaculture growth, S. latissima blades can grow between 60-140 cm depending on the water depth, planting time, and nutrient availability (Handå et al., 2013). Oysters, however, are the most widely aquacultured species in coastal areas of the U.S (NMFS, 2018). The Eastern oyster (Crassostrea virginica) mostly grows during the summer months where water temperatures are above 15 °C and is in a state of relative dormancy in the winter (Dame, 1972). Therefore, it has been suggested that S. latissima could complement farmed oysters because of the differences in growing season, thus providing an additional source of income without interfering with oyster production. This new industry, however, would benefit from production estimates that would enable farmers to decide whether it is possible or lucrative to farm kelp.
Bioenergetics allows for study of factors influencing an individual species' growth through assessment of energy relationships and transformations. Thus, bioenergetics models can support the creation of production estimates. Growth of S.
latissima has been studied in both field and laboratory settings in response to environmental conditions. For example, the impacts of irradiance, temperature, and nutrient concentration on S. latissima growth were highlighted through a transplant experiment with kelps at different depths (Boden, 1979). In a simple predictive model created for S. latissima, growth was assumed to have a linear relationship with dissolved inorganic nitrogen concentration (with a temperature correction; Petrell et al., 1993). This model required an assumption that nitrogen dynamics are always limiting growth, thus ignoring the potential influence of irradiance. Light-harvesting characteristics of S. latissima have been shown to vary in populations from different ambient light regimes (Gerard, 1988), so irradiance is an important factor on this species' dynamics. Dynamic energy budget (DEB) theory provides a framework to examine the interactive effects of environmental nutrient concentrations and irradiance on an organism through parallel systems of nitrogen and carbon dynamics (Kooijman, 2010). DEB theory is a formal theory of metabolism, which is used to model the flow of mass and energy through an organism from uptake to usage for growth, maintenance, reproduction, or excretion. DEB theory provides a sound mechanistic basis for understanding an organism's energetics.
Modeling autotrophs is a relatively new direction for DEB theory. The standard DEB model, which has been applied to many animal species, is appropriate for an animal that does not change shape through development and eats only one food source with a constant chemical composition. The standard DEB model has one reserve and one structure; structure and reserve are the central pools of matter modeled as state variables within an organism. For the application of DEB theory to autotrophs, such as S. latissima, multiple reserves are necessary to accurately model matter and energy dynamics because nutrient uptake can occur separately from photosynthesis (Kooijman, 2010). Autotroph DEB models have been constructed for microalgae (Lorena et al., 2010, Livanou et al., 2019, phytoplankton-zooplankton interactions (Poggiale et al., 2010), calcification of a coccolithophore (Muller and Nisbet, 2014), and the macroalga Ulva lactuca (Lavaud et al., under review). Other modelers have used dynamic bioenergetic models that borrow some concepts from DEB theory but selected a greater degree of simplification to model S. latissima (Broch and Slagstad, 2012) and the interaction between heterotrophic coral and autotrophic Symbiodinium spp. (Cunning et al., 2017). The model of Norwegian S. latissima was created as a tool for optimizing aquaculture production (Broch and Slagstad, 2012). They chose a simpler base structure than Lorena et al. (2010), but this simplification does not increase parsimony (i.e. reduce the number of model parameters) because of their use of correction functions where a DEB model would not need them. For example, a simple functional response is used to control the effect of frond size on gross growth rate because smaller kelps grow faster than larger kelps (Broch and Slagstad, 2012).
This would be redundant in a DEB model with maintenance because the maintenance cost is volume specific (Lorena et al., 2010).
Our primary objective is to develop a bioenergetics model for S. latissima growth. Specifically, we aim to calibrate the general structure of a macroalgae DEB model to our field data from Rhode Island (U.S.). The resulting model allows for growth predictions based on inputs of environmental conditions and has the potential to support the sustainable aquaculture industry, particularly with regard to site selection.

Dynamic Energy Budget Model Structure
The core structure of the S. latissima model tracked the uptake of carbon and nitrogen, their assimilation into reserves and allocation to growth, maintenance, or excretion ( Figure 1). The variables that depict the state of the model were nitrogen reserve density, carbon reserve density, and structure. There were three differential equations to represent the change in the state variables over time. The differential equation solver used to run the model in R (R Core Team, 2019) was the package deSolve . The model initial conditions (Table S1a) and molecular weights for the structure and the reserves (  (Schiener et al., 2015).
Therefore, we modeled the carbon reserves with a chemical composition similar to laminarin, mannitol, and glucose setting a constant ratio of one carbon to two hydrogen to one oxygen (1:2:1, C:H:O).
For our S. latissima DEB model, we considered two reserve pools: carbon and nitrogen (nitrate and ammonium, collectively); other potential nutrients such as phosphorous or potassium were dismissed. We feel this is a valid assumption, especially in regions where nitrogen is not abundant year-round and nitrogen availability is what drives accelerated growth in winter and early spring (Gagné et al., 1982). A previous dynamic model for Norwegian S. latissima growth also chose to use only a carbon and a nitrogen reserve (Broch and Slagstad, 2012). Adding further reserves to the model would increase complexity by increasing the number of state variables and parameters with potentially little to no accompanying increase in accuracy. Our model also assumed that the energy carbon concentrating mechanisms used to convert bicarbonate into carbon dioxide do not have a significant impact on the overall energetics. The majority of the dissolved inorganic carbon in the ocean is in the form of bicarbonate (Raven et al., 2005); S. latissima and other algae use carbonic anhydrase in carbon concentrating mechanisms to assimilate bicarbonate and convert it into carbon dioxide (Axelsson et al., 2000). In other words, we assumed that assimilating carbon dioxide directly was identical to assimilating carbon dioxide that was formed extracellularly from bicarbonate through a carbon concentrating mechanism.
Another assumption of our DEB model was that S. latissima can be considered as a V1-morph. In DEB theory, V1-morphs are organisms whose surface area is proportional to volume, which simplifies the dynamics of the model (Kooijman, 2010). Saccharina latissima grows as a sheet in both length and width directions at the meristematic blade region near the stipe (Sjøtun, 1993). We were assuming variation in blade thickness over an individual blade and through time (Vettori and Nikora, 2017) would not have a substantial enough impact on this surface area to volume ratio to preclude the V-1 morph assumption. Assuming that surface area was proportional to volume is reasonable at sites where there are not major differences in water velocity.
Drag from water speed has been found to change blade morphology; more sheltered environments have led to wider blades with ruffled edges and areas with strong currents have caused blades to be more strap-like (Buck and Buchholz, 2005). Even with this plastic morphology, the change in blade height was still small enough to not significantly impact surface areas in regions with similar water speeds.
Other assumptions for this model were grounded in reducing the time of year that the model can be logically applied. Energy was not used for reproduction or maturity in our model for S. latissima as it would be in a standard DEB model, which is a simplification that allows for a more parsimonious model. Saccharina latissima peaks in reproduction twice a year, June and October, but can have reproductive tissue year-round (Lee and Brinkhuis, 1986). There is evidence it produces inhibitors that minimize the formation of reproductive tissue during the rapid growth phase Lüning, 1999, Lüning et al., 2000). Reproductive development happens, but only for a small subset of blades by the time the aquaculture harvest occurs in Spring towards the end of first period of rapid growth. Apical frond loss in kelps is correlated with temperature stress and wave action (Krumhansl et al., 2014), mechanical stress of biofouling (Brown et al., 1997), and overall blade length (Sjøtun, 1993). Our model does not incorporate a correction for apical frond loss because we were focusing on the aquaculture season and the exact mechanism for this loss remains very context-specific in the literature. Furthermore, aquaculture farmers generally harvest kelp before biofouling begins, which maximizes harvestable blade length. We also assumed photoinhibition does not occur, which is similarly reasonable for the aquaculture season but not the summer. Photoinhibition does occur in S.
latissima especially when high light conditions are combined with high temperature conditions (Heinrich et al., 2012). The aquaculture season of kelp is placed to maximize growth while minimizing loss or degradation of tissues due to various stresses, and a tissue loss function would be necessary to accurately model wild kelp year-round.

Empirical Data
Saccharina latissima was grown at four oyster farm sites from fall to spring in both 2017-2018 and 2018-2019 (see Table S3 for exact plant and harvest dates). Kelp seed was raised in aquaria from harvested local reproductive kelp tissue collected at

Environmental Forcing Functions
The S. latissima model was forced with temperature, irradiance, dissolved inorganic carbon (DIC) concentration, and nitrate and ammonium concentration data on an hourly time step. Temperature recorded at fifteen-minute intervals were averaged on an hourly basis and was identical at each kelp line at a specific site and year in the model because we only had a logger on one of the longlines at each site and assumed the conditions were the same at the other. Due to difficulties with biofouling on irradiance loggers, we used radiative forcing from the North American Regional Reanalysis (Mesinger et al., 2006) to estimate photosynthetically active radiation (PAR; Appendix S4). We used linear interpolation to create an hourly forcing from source data from every three hours ( Figure 3). All sites have the same base irradiance forcing in one year using this method. DIC concentration data were not collected in this study, so this forcing is estimated from other sources. The Pt. Judith Pond sites were held at a constant DIC value based on U.S. Environmental Protection Agency data from Ninigret Pond (J. Grear, unpublished data). The Narragansett Bay sites were held at a constant DIC value based on data from Brenton Point (Segarra, 2002). We used linear interpolation to estimate hourly nitrate and ammonium concentrations.

Model Calibration
The parameters for our model were set using a combination of manual fitting to literature data, values from previous autotroph DEB models, and field data from this study (Table 1). Literature data on basic physiological rates were compiled to simultaneously calibrate parameters using the model equations (Table 3). Information about the locations where these studies were conducted was also included because there are multiple ecotypes of S. latissima (Gerard, 1988), which may influence their physiological response. Due to a lack of local information on certain aspects of kelp growth, this model was calibrated with data across multiple ecotypes of kelp. The Arrhenius relationship parameters were calculated using a least squared non-linear regression on the compiled literature rates which were standardized using the reference temperature. The nitrate and ammonium uptake parameters, maximum volume specific nitrogen assimilation and half-saturation concentration for NO3and NH4 + uptake, were calibrated using nitrate uptake data from Espinoza and Chapman (1983). Photosynthesis parameters, photosynthetic unit (PSU) density, binding probability of a photon to a free light synthesizing unit, and dissociation rate of releasing ATP and NADPH + , were calibrated using oxygen production data from Johansson and Snoeijs (2002). Root mean square error (RMSE) was used as a measure of spread in the residuals for assessing the quality of each parameter calibration.

Sensitivity Analyses
To determine how each DEB parameter influenced results, we analyzed the local sensitivity of C reserve density, N reserve density, and structural mass to model parameters using an L1 summary value of sensitivity from the R package FME . L1 is a summation of the absolute value of all the elements in a sensitivity matrix divided by the number of elements. Each element of the sensitivity matrix is calculated by multiplying the change in the output variable over the change in a parameter to the scaling of that same parameter over the scaling of that same output variable. The scaling here is usually equal to the value of the variable or parameter . Parameters have unequal impacts on model outcomes, and this method allowed us to determine what parameters are having the largest effects on , , and .

Model Calibration: Literature Data
The Arrhenius relationship fit to the compiled literature data (Table 3) reflected maximum physiological rates at temperatures around 13 °C ( Figure 4). The lower boundary of the Arrhenius relationship was 0 °C, and the upper boundary was 13.39 °C. The adjusted R-squared for this relationship was 0.55 (p-value = 2.74e-11).
Using the nitrate uptake data from Espinoza and Chapman's (1983) site that was seasonally depleted of nitrate between April and November provided estimates of maximum volume specific ammonium and nitrate assimilation of 2.7 *10 -4 (mol N mol MV -1 h -1 ) and a half-saturation concentration of 2.667 *10 -6 (mol NO3and NH4 + L -1 ) ( Figure 5). The fit for the data collected at 18 °C was slightly better at a RMSE of 9.73e-07 (mol N gDW -1 h -1 ) than the 9 °C data at 1.43e-06 (mol N gDW -1 h -1 ).
For the oxygen production data (Johansson and Snoeijs, 2002) used to calibrate the photosynthesis parameters, the values that had the lowest error around the data were a ()* photosynthetic unit (PSU) density of 0.05 (mol PSU mol MV -1 ), a ̇binding probability of a photon to a free light synthesizing unit of 2.8 *10 -6 (unitless), and a ̇dissociation rate of releasing ATP and NADPH + of 0.28 (mol NADPH mol PSU -1 h -1 ; Figure 6). The resulting RMSE was 0.000457 (g O2 g DW -1 h -1 ). The maximum oxygen production rate of the model was around 0.00495 g O2 g DW -1 h -1 ( Figure 6).
Appropriate literature data for calibrating carbon dioxide uptake were not available. Air-based carbon dioxide uptake data for S. latissima (Ní Longphuirt et al., 2013) were examined but ultimately rejected due to likely dissimilarity to submerged carbon dioxide uptake. Other parameters in DEB models such as the reserve turnover rates could not be compared to measurable physiological data due to the internal nature of the processes without easily interpretable signals, so these parameters were set at values similar to those in Lorena et al. (2010) and Lavaud et al. (under review). Since both irradiance and DIC concentration impact carbon assimilation, plots of these initial assimilation fluxes revealed which one of these factors limits carbon assimilation. The specific assimilation rate of carbon mirrored the shape of the specific relaxation rate ( Figure S5c and S5d). The shape of the specific relaxation rate ( Figure   S5c) reflected the influence of the temperature correction ( Figure 7) rather than that of the magnitude of the irradiance forcing ( Figure 3).

Sensitivity Analyses
The parameters with the largest effects (>3000 L1 summary value of sensitivity functions) on the state variables were C reserve density, N reserve density, and structural mass were reference temperature, upper boundary of temperature tolerance, lower boundary of temperature tolerance, Arrhenius temperature outside TH, Arrhenius temperature outside TL, yield factor of C reserve to structure, fraction of rejection flux incorporated back in ireserve, yield factor of C reserve to NADPH, and yield factor of C reserve to CO2 (Figure 11). The parameters with the moderate effects (1000-2000) on the state variables were ̇ max. volume specific carbon assimilation, and ̇ dissociation rate of releasing ATP and NADPH + . The parameters with small but non-zero effects on the state variables were yield factor of N reserve to structure, ̇ carbon reserve turnover rate, and photosynthetic unit (PSU) density ( Figure 11).

Discussion
Further mariculture development represents a key role in expanding U.S.
sustainable food production, and macroalgae can provide high returns if the proper growth conditions exist. Thus, understanding the growth dynamics of S. latissima can provide the industry with a powerful predictive tool for estimating production. Our model is the first attempt to apply Dynamic Energy Budget (DEB) theory to a macroalga of the order Laminariales. The model tended to underestimate initial growth of S. latissima due to a low modeled specific relaxation rate or for sites where N was limiting an inaccurate N forcing. The model was good, however, at predicting final growth rates ( Figure 10). Our process-based model also allowed us to better understand local growth limitations as they relate to the behavior of the model.

Growth Limitation
We hypothesize that the early season model underestimation of S. latissima length growth results from limited modeled carbon assimilation due to a low modeled specific relaxation rate (i.e. the light-dependent reactions of photosynthesis) for those model runs where C dynamics are the limiting factor. Seasonal temperature change and diurnal irradiance oscillation appear to control C assimilation via the specific relaxation rate rather than a seasonal trend of irradiance magnitude change. There is some evidence that winter limitation of C dynamics occurs in the field; S. latissima individuals older than a year were shown to have a decrease in blade C content midwinter suggesting consumption of stored carbohydrates (Sjøtun, 1993). New sporophytes would not have this carbohydrate pool to draw upon. This decrease in C content suggests that C dynamics may be limiting S. latissima growth, but limitation was not directly examined by Sjøtun (1993).
For those model runs where C dynamics are the limiting factor, the early season (planting to end of March) model underestimation of field growth may reflect that the temperature correction has too powerful of an impact on the specific relaxation rate in comparison to irradiance. The seasonal trend of the specific relaxation rate reflects the seasonality of temperature; we suspect that the assumption of temperature impacting ̇dissociation rate of releasing ATP and NADPH + may cause temperature to have an outsized impact in comparison to irradiance magnitude.
Ocean temperature trends trail behind irradiance change, and Narragansett Bay is no exception to this pattern (Brady-Campbell et al., 1984). Saccharina latissima's early season growth could be driven by this early season increase in irradiance rather than water temperature change. This possibility is further reason to verify an organism's Arrhenius relationship regionally if there is potential for different degrees of temperature adaptation. More data is necessary to confirm if the S. latissima model underestimates winter kelp growth due to the impact of temperature on modeled photosynthesis.
Other than an increase in irradiance, daylength is a linked variable that could impact seasonal growth patterns. Broch and Slagstad's (2012) S. latissima model used the rate of change of day length in a photoperiodic effect function to create seasonality in their growth prediction. This is based on the hypothesis that S. latissima is a "seasonal anticipator" with endogenous circadian rhythms (Kain, 1989). Seasonal anticipators are posited to grow strategically in response to a trigger as opposed to simply responding to environmental conditions. Other kelps, Laminaria hyperborea and Laminaria digitata, have been shown to have free-running seasonal growth patterns which suggests control by endogenous circadian rhythms (Schaffelke and Lüning, 1994 (Kremer and Markham, 1979). The linkage between the lightdependent and light-independent reactions is modeled as an immediate transference. In other words, when there is no irradiance input to the S. latissima model, the assimilation of carbohydrates to the carbon reserve is zero. This is not reflective of the lag which allows for the formation of carbohydrates in the dark, but we would argue that adding this layer of physiological accuracy would reduce model efficiency without enough increase in the predictive capacity.
Where N dynamics limit growth below the level observed in the field, the heavily interpolated N forcing is likely not reflective of actual conditions. The N forcing effects our ability to fit the overall model tightly to all runs simultaneously due to overestimation of blade length in some runs with an overinflated N forcing and underestimation in others.

Sensitivity Analysis
Identifying parameters that a model is highly sensitive to gives us a clearer idea of what parameters are likely to shift the overall fit of a model . The high sensitivity of the state variables to the temperature related parameters is a logical outcome of the central role of temperature in DEB theory.
Since the temperature correction is applied to such a large number of rates in the organism, the high sensitivity to these values is reasonable. The sensitivity of the model to the temperature parameters is an argument for caution in regional calibration  (Gerard, 1988). For instance, S. latissima individuals from New York have been shown to have a different physiological response to temperature in a lab setting than individuals from Maine (Gerard, 1988). Narragansett Bay, RI is located towards the southern boundary of where S. latissima can survive (Taylor, 1972). The existence of multiple ecotypes of this species suggests that parameters may require regional adjustment, particularly in the Arctic. Our model assumptions about surface area's proportionality to volume impedes prediction of blade shape plasticity, which is a characteristic of S. latissima related to drag (Buck and Buchholz, 2005). Since the blade thickness and amount of blade ruffling impact the relationship between surface area and volume, this complicates the use of the model in areas with either significantly different amounts of drag from water motion or highly variable water motion.
Further research on the mechanisms for frond loss, blade plasticity, and regional parameter information have the potential to improve this DEB model. A clearer physiological cause for apical frond loss would allow this to be included in mechanistic models in a more meaningful way than a correction function setting erosion based on one correlated variable like length or age. Similarly, determining a mechanism for how blade type changes in response to variable water speeds would provide a clearer picture of overall growth dynamics. Other data that would be useful for understanding S. latissima physiology and increasing the accuracy of predictive modeling are underwater carbon dioxide uptake data in response to variable irradiance and more regionally appropriate oxygen production data in response to variable irradiance.
Our S. latissima model is a first step towards estimating kelp aquaculture production in the U.S. Norwegian S. latissima research highlights the potential of dynamic modeling to both estimate the effectiveness of S. latissima as a nutrient assimilator in integrated multi-trophic aquaculture (IMTA) systems , Fossberg et al., 2018 and estimate aquaculture production at large scales using a coupled biogeochemical-hydrodynamic-kelp model (Broch et al., 2019). In future work, our S. latissima DEB model could be coupled with a DEB model for C.
virginica (Filgueira et al., 2014;Lavaud et al., 2017) and the Regional Ocean Modeling System (ROMS) with a Carbon Silicate Nitrogen Ecosystem (CoSiNE) model (Chai et al, 2009) to predict growth potential at sites. ROMS is a widely-used primitive equations ocean model and the CoSINE model integrates biogeochemical processes (Chai et al, 2009 Table   Table 2. Model equations with environmental conditions: T = temperature (K), I = irradiance (μE m -2 h -1 ), DIC = dissolved inorganic carbon (mol DIC L -1 ), and N = nitrate and ammonium concentration (NO3and NH4 + L -1 ).       . The Arrhenius relationship for S. latissima was estimated using multiple growth and photosynthesis datasets: 1) Bolton and Lüning (1982) (squares; black for kelp from France, dark grey for Norway, light grey for Germany, white for the UK), 2) Fortes and Lüning (1980) (diamond), 3) Davison and Davison (1987) (asterisk), and 4) Davison (1987) (circles; black for sporophyte rearing temp 0°C, dark grey for 5°C, dark grey with black border for 10°C, light grey for 15°C, and white for 20°C). The adjusted R-squared statistic for the fit of the curve to the data points is 0.551 (p-value = 2.74e-11).  Espinoza and Chapman (1983). Black depicts the uptake at 9°C and grey illustrates the uptake at 18°C. The RMSE for the model calibrated to the 9°C data is 1.43e-06 and 9.73e-07 (both mol N gDW -1 h -1 ) for the model calibrated to the data collected at 18°C. . Modeled (lines) and observed (dots) oxygen production from Johansson and Snoeijs (2002). The RMSE for the fit of this model curve to this data is 0.000457 (g O2 g DW -1 h -1 ).

SUPPORTING INFORMATION
Appendix S1: Model Initial Conditions and Molecular Weights Initial reserve densities were set based on the knowledge that the carbon reserve would occupy a larger amount of mass in comparison to the nitrogen reserve at the starting time of the study. Initial blade lengths and biomass were set together using the Gevaert et al. (2001) power relationship to set a length similar to that of the sporophytes on the day they were planted out on the farms. Initial mass of structure was set to calculate the chosen initial blade biomass. The molecular weights used in this biomass calculation were calculated from assumed C:H:O:N chemical compositions and the associated atomic masses of these elements. The nitrogen reserve's chemical composition is ammonia, the carbon reserve's chemical composition is glucose, laminarin, and mannitol, and the structure's chemical composition is the average of mannitol and laminarin. Because structure is also composed of a significant amount of proteins, we set the structure chemical index for N to 0.04.

Appendix S2: Model Equation Descriptions
A temperature correction (CT) is applied to all modelled rates expect for the photon binding rate (̇-) (Lorena et al., 2010, Lavaud et al., under review). Arrhenius relationships are based on the idea that metabolic rates of a species are impacted the same amount by temperature change (Kooijman, 2010). This correction is an extended Arrhenius correction factor; the extension supports a non-linear response to temperature change: With TA the Arrhenius temperature, T0 the reference temperature, T is the field or input temperature, TL and TH the lower and higher bounds of temperature tolerance, and TAL and TAH the Arrhenius temperatures for the lower and higher temperatures outside the optimal range (all in K).
In this model, rates (variables with dotted letters) are expressed per C-mol of structure. The specific assimilation rate of nitrate and ammonium (̇~• € in mol N mol MV -1 h -1 ) is calculated through a Michaelis-Menten relationship:

̇=̇ * * [ ]/([ ] + )
With ̇~• €• maximum volume specific nitrogen assimilation (mol N mol MV -1 h -1 ), [N] the concentration of ammonium and nitrate in the environment (NO3and NH4 + L -enzyme kinetics with a major difference being the assumption of low substrateenzyme dissociation rates (Kooijman, 2010). SUs allow for dynamic substrate limitation and flexible computation. For carbon dynamics in particular three linked SUs model the different parts of photosynthesis. One of these SUs facilitates the uptake of carbon dioxide. The specific CO2 uptake rate (‚ ƒ " in mol CO2 mol MV -1 h -1 ) is calculated through a Michaelis-Menten relationship similarly to ̇~• € in this model: With ̇‚ ƒ " • maximum volume specific CO2 uptake rate (mol CO2 mol MV -1 h -1 ), [DIC] the concentration of dissolved inorganic carbon (mol DIC L -1 ), and KC half-saturation concentration for CO2 uptake (mol DIC L -1 ). Bicarbonate is converted to carbon dioxide by carbonic anhydrase, so this is why the forcing is DIC and what is being taken in is CO2.
Another of the photosynthesis SUs depicts the light-dependent reactions of photosynthesis. In this equation, photons are bound to photosystems and ATP and NADPH are produced and released at set rates. The specific relaxation rate (-in mol NADPH mol MV -1 h -1 ) is given by: ̇= * * + * * With ()* photosynthetic unit (PSU) density (mol PSU mol MV -1 ), I forcing function irradiance (μE m -2 h -1 ), ̇photon binding rate (unitless), and ̇dissociation rate of releasing ATP and NADPH + (mol NADPH mol PSU -1 h -1 ).
The final SU of the carbon dynamics integrates the inputs of the previous two SUs and assimilates carbon to its reserve. This part of the photosynthesis model represents the Calvin-Benson cycle during which carbohydrates are formed from carbon dioxide. The specific assimilation rate of C (̇ † € in mol C mol MV -1 h -1 ) is calculated through the parallel processing of NADPH and CO2: With ̇~ † €• maximum volume specific carbon assimilation (mol C mol MV -1 h -1 ), ‚ƒ " ‚ yield factor of C reserve to CO2 (mol CO2 mol EC -1 ), and -‚ yield factor of C reserve to NADPH (mol NADPH mol EC -1 ).
Oxygen is produced as a byproduct of the reduction of NADP+, and one mole of O2 is produced for every four moles of NADPH produced. Calculating the oxygen production rate (ƒ " in g O2 g biomass -1 h -1 ) is a useful equation in this model to provide comparisons to literature data to calibrate the photosynthesis parameters: ̇=̇ * * * With MV structural mass (mol MV), ƒ " molar weight of O2 (g O2 mol O2 -1 ), and B is modeled kelp blade biomass (g).
Reserves are catabolized to send mass on to be used in growth or maintenance.
The specific catabolic flux of the N or C reserves (~Š ‚ in mol EN or EC molMV -1 h -1 ) follows first order kinetics, which means that as the carbon or nitrogen reserve density (~Š in mol EN or EC mol MV -1 ) increases the faster catabolism will occur: = ]^ * ) −̇` With ~OE carbon or nitrogen reserve turnover rate (h -1 ) and ̇ net specific growth rate (h -1 ):

̇=
̇ is modified at each timestep using a function modified from the DEBtool package's function called sgr2.m in the alga section (https://www.bio.vu.nl/thb/deb/deblab/debtool/DEBtool_M/manual/index.html). The method that underlies this function is the Newton-Raphson method with continuation.
Maintenance takes priority over growth for catabolized resources in DEB theory. The specific flux for metabolism from the N or C reserves (̇~Š Ž Š in mol EN or EC mol Mv -1 h -1 ) is set by the temperature corrected volume-specific maintenance cost paid by N or C reserve (~Š Ž in mol EN or EC mol Mv -1 h -1 ) if the catabolic flux is enough to meet those costs: The specific growth fluxes from both N and C reserves (~Š • in mol EN or EC mol Mv -1 h -1 ) sent to the growth SU contain any mass left after maintenance costs are filled:

= −
If maintenance requirements not met by the catabolic fluxes from the reserves, maintenance costs will be met by drawing from structure. The specific maintenance flux from structure ( • Ž in mol N and C mol Mv -1 h -1 ) is given by: With ~Š• yield factor of N reserve or C reserve to structure (mol EN or EC mol Mv -1 ).
The growth SU handles growth fluxes from both reserves complementarily, so the specific gross growth rate ( •• in h -1 ) is given by: Either the amount of carbon or nitrogen available will limit the growth SU, so excess will be available for one of the compounds. The model can reject this excess material as binding sites are filled and return it to the respective reserve. This rejected specific C or N flux from growth SU (~Š ' in mol Ei mol MV -1 h -1 ) is the difference between available growth flux and what is actually used for growth:

= − *
The dynamics of the state variables are the main differential equations of the model. The dynamics of the N or C reserve densities balance the inputs of assimilation and returning rejected flux with the outputs of catabolism and dilution by growth: With ~Š fraction of rejection flux incorporated in i-reserve (unitless).
The dynamics of structural mass are controlled by the net specific growth rate:

=̇ *
Modeled kelp blade biomass (dry weight) is calculated by summing the mass of the reserves and the structure in the model: With molar weight of structure (g mol -1 ), molar weight of C reserve (g C mol C -1 ), and molar weight of N reserve (g N mol N -1 ).
Kelp blade length (cm) is calculated using an allometric relationship between length and dry weight from Gevaert et al. (2001): = .