INCORPORATING TEMPERATURE-DEPENDENT FISH BIOENERGETICS INTO A NARRAGANSETT BAY FOOD WEB MODEL

Food web models capture shifting species interactions, making them useful tools for exploring community responses to perturbations. The inclusion of environmental drivers, such as temperature, can improve model predictions as energy demands of an organism can be temperature-specific. While Ecopath with Ecosim (EwE) and the recent R implementation of this software, Rpath, have included some thermal responses in past work, models have yet to include temperature-dependent energetic demands and metabolic costs. Our work demonstrates the inclusion of temperaturedependent bioenergetics into an Rpath food web model using the case study of a warming estuary: Narragansett Bay (RI, U.S.). Thermal response parameters from literature were used to construct Kitchell curves describing temperature-dependent consumption and modified Arrhenius curves describing temperature-dependent respiration. Surface water temperature time series from 1994 to 2054 for high and low warming scenarios were created from observed temperatures and projections from the Coupled Model Intercomparison Project (CMIP6) multi-model ensemble. The integration of temperature-dependent fish bioenergetics resulted in lower projected biomasses compared to the base version of the model without environmental forcing, reflecting the impact of increased energetic demands. The differences in the modelpredicted biomasses highlight the importance of accounting for thermal effects on marine species in ecosystem models, which will become increasingly important as ocean temperatures continue to rise in Narragansett Bay and worldwide.

. Biomass, the fraction of energy towards respiration (ActiveRespFrac), and the  reflecting the impact of increased energetic demands. The differences in the modelpredicted biomasses highlight the importance of accounting for thermal effects on marine species in ecosystem models, which will become increasingly important as ocean temperatures continue to rise in Narragansett Bay and worldwide.

Introduction
Climate change represents a major continuous perturbation to the Northeast U.S. Continental Shelf ecosystem, affecting species distribution, migration phenology, and physiological processes Kleisner et al., 2017;Langan et al., 2021;Pershing et al., 2015). Rises in sea temperature have been acutely observed in Narragansett Bay, Rhode Island (U.S.), where surface water temperatures have risen approximately 1.5°C since 1950 (Fulweiler et al., 2015). Climate-driven models can help elucidate how these changes in temperature are affecting ecosystem dynamics (Brander, 2015). Food web models, such as those created using the Ecopath with Ecosim (EwE) software, are becoming increasingly popular to study how ecosystems respond to changes in fisheries harvest, species interactions, and external drivers (Buchheister et al., 2017;Colléter et al., 2015;Villasante et al., 2016). While some studies have included environmental components (Bentley et al., 2017;Corrales et al., 2017;Serpetti et al., 2017), there are thermal influences on marine populations, such as direct impacts on species bioenergetic demands, that have not yet been incorporated into EwE models.
As ectotherms, fish rely on the environment to regulate their body temperature, with ambient temperatures ultimately influencing physiological rates . Consequently, bioenergetics and physiological responses are often used as the basis for mechanistically-driven models (Jørgensen et al. 2016). The basics of temperaturedependent bioenergetics have been recognized for decades (Brett, 1971;, and while species-specific rates are often unknown, it is nonetheless important to begin incorporating bioenergetic relationships into multispecies ecosystem models for increased realism of how species dynamics are affected by warming waters. Understanding community responses to climate change, rather than examining singlespecies responses in isolation, will provide more insight as to how environmental stressors will affect marine ecosystems (Nagelkerken and Munday 2016).
As temperatures increase from the cooler portion of a species' thermal tolerance, it can become more energetically expensive for ectotherms to maintain base metabolic demands . The thermal response of metabolism is frequently described with an Arrhenius equation (Eq. (1)) if thermodynamic relationships are considered the dominant drivers (Brown et al., 2004;Gillooly et al., 2001;Schulte, 2015). The Arrhenius equation, often used in fitting laboratory data or in models incorporating metabolic ecology (Blanchard et al., 2012;Clarke and Johnston, 1999;Dahlke et al., 2020;Neubauer and Andersen, 2019), calculates the rate of reaction (k) as a function of a constant (A), the activation energy (Ea), the universal gas constant (R), and the temperature (T) in degrees kelvin.
Temperature-dependent energetic costs can modify species interactions, primarily through the adjustment of consumption rates as predators alter their intake to maintain their energy balance (Johansen et al., 2015). In relation to temperature, consumption increases as waters warm, reaches a maximum at some optimum temperature, and then ingestion sharply decreases as the maximum tolerated temperature is approached (Fogarty and Collie, 2020;. Increasing consumption with increasing temperatures has been documented for a variety of fishes including bluefish (Pomatomus saltatrix), brook trout (Salvelinus fontinalis), and anemonefish (Amphiprion melanopus) Nowicki et al., 2012;Olla et al., 1985;Ries and Perry, 1995). Increased metabolic demand in warmer waters likely accounts for the increasing portion of the curve, while stress responses and behavioral shifts are thought to drive the reduced consumption rates frequently seen near species' thermal maxima (Brett, 1971;Jobling, 1997;Johansen et al., 2014;Nowicki et al., 2012). The exact shape can vary due to the many factors affecting consumption, such as locomotion, hormone regulators, detection, and successful prey capture, but numerous experimental studies have documented this general thermal response (Jobling, 1997;Volkoff and Rønnestad, 2020).
Higher energetic demands can adversely affect production, or the surplus energy available for optional processes such as growth and reproduction Neubauer and Andersen, 2019). Temperature can adjust the efficiency with which organisms transform food energy into growth (Lemoine and Burkepile, 2012).
Given that the foundational bioenergetic relationships apply to many marine species operating in their preferred thermal range (Deslauriers et al., 2017;Sibly et al., 2012), the integration of temperature-dependent bioenergetics in ecosystem models can provide more realistic predictions of how climate change will impact ecosystem production .
The EwE software has been used to model over 400 marine and aquatic ecosystems (Colléter et al., 2015). Ecopath creates static, mass-balance food web models that give a snapshot of the energy flow of an ecosystem (Polovina, 1984).
Ecopath is governed by master equations for consumption and production and requires data inputs of biomass (B), production to biomass ratio (P/B), consumption to biomass ratio (Q/B), a diet matrix, and fishing information (Christensen and Pauly, 1992).
Ecosim expands this snapshot to create time-dynamic projections of biomass (Coll et al., 2009). Predation is modeled with foraging arena theory in which prey are, at times, vulnerable to predation and invulnerable at others (Ahrens et al., 2012;Walters and Christensen, 2007). The vulnerabilities, or the parameters specifying the rate of exchange between vulnerable and invulnerable states, are estimated with a fitting procedure to minimize the sum of squares between projected and observed biomasses . External forcing functions of changing inputs, such as primary production or fishing effort, can also be used to drive the dynamic models into the future (Christensen and Walters, 2004). with temperature-dependent consumption used temperature of occurrence to estimate thermal preference, and modified consumption to restrict foraging capacity for each species or functional group. However, no EwE model has so far incorporated the other major energetic impact of temperature: changing energy demands and metabolic costs.
One of the master equations of Ecopath specifies that consumption in units of biomass consumed is the sum of production, respiration, and unassimilated food (C = P + R + U ; Christensen et al. 2005). EwE aggregates metabolic costs into the respiration term which represents the biomass lost to the ecosystem (i.e. energy consumed and assimilated by a group but is not transformed into production). Thermal adjustments to the respiration term are not in the current EwE software capacity. However, EwE-type food web models can now be created and investigated in the flexible R software using the 'Rpath' package (Lucey et al., 2020). The Rpath package creates an opportunity for further exploration of the ecosystem responds to shifting environmental conditions and the inclusion of temperature-dependent bioenergetics.
The goal of this study is to demonstrate how temperature-dependent bioenergetics can be incorporated in a food web model using Rpath, the flexible R implementation of the Ecopath with Ecosim modeling framework. We apply this method to better understand climate change impacts using a rapidly warming Narragansett Bay as a case study. We hope to expand the functionality of Rpath and illustrate the scales at which temperature-dependent consumption and respiration can amplify to alter food web model biomass outputs.  Table A1.2) that compose the fish functional groups.

Temperature
Surface water temperatures were taken from the University of Rhode Island Graduate School of Oceanography (URI GSO) weekly fish trawl (Collie et al., 2008). through range shift, ecological feedbacks, and evolution (Peck, 2011;Sibly et al., 2012). The final time series for the high and low warming scenarios were created by averaging the output of each of the six delta corrected surface temperature projections for each year after 2018.

Temperature-dependent consumption
The food consumption thermal response for each of the three fish functional groups was described using the Kitchell equation Kitchell et al., 1977), which is routinely used to characterize temperature-dependent consumption in bioenergetic models (Hansson et al., 1996;. The Kitchell equation, shown in Eq. (2), uses straightforward input parameters of the thermal maximum (Tmax; the temperature above which consumption is zero), temperature of optimum consumption (ToptC), and the Q10 of consumption (referring to the rate of change for a process as the temperature increases 10°C; , Fogarty & Collie 2020 to estimate the proportion of maximum consumption that occurs at a given temperature (rc). We ran a sensitivity test to examine the effect of input community composition on the consumption thermal response curves. Three Kitchell curves per functional group, weighted by different species biomasses, were used to create consumption modifier time series; 1) the 1994-1998 averaged biomasses as described earlier in the methods, 2) the single year's observed biomass from the time series data that resulted in the strongest warm-skewed curve, and 3) biomasses including more southern, warm-water species to represent a future curve as new species enter the Bay. Further information on the creation of the curves with southern species can be found in Supplement ( 3)). The scaled Kitchell curve, referred to as the relative consumption curve (RelC), had a maximum consumption modifier greater than one at ToptC.
The base model with temperature-dependent consumption will be referred to as the 'consumption' version of the model. Temperature-dependent consumption was forced differently than previous EwE models because the forcing functions between Ecosim and Rsim differ. In this study, the consumption modifiers, as calculated from the relative consumption curve and temperature time series, were applied with the 'ForcedSearch' feature (formerly called ForcedPred). This forcing option modifies the effective predator biomass which is then inputted into the main consumption calculations (Lucey et al. 2020 equations 19 & 20). The calculations of ForcedSearch function apply the thermal response before the full consumption calculations, so that the foraging arena and interspecies interactions mediate the temperature effect on physiological processes, as has been suggested by others (Neubauer and Andersen, 2019). A time series of consumption modifiers by year then forced the temperaturedependent model versions using the 'adjust.forcing' function on the ForcedSearch parameter.

Temperature-dependent respiration
The second bioenergetic response built into Rpath was temperature-dependent respiration. Respiration was treated similarly to standard metabolism, though respiration includes other energetic costs including standard dynamic action (SDA), the metabolic cost to digest food, or activity costs. Our framing of metabolism as the amount of energy used to maintain function was more directly applicable to the EwE mass-balance setup and respiration term than other metrics such as aerobic scope. The original Innes-Gold et al. (2020) EwE model was a parsimonious model of the ecosystem in which species age and size were generally not included. Therefore, our functional form for temperature-dependent respiration was simplified because it did not include a body mass effect on respiration (Sibly et al., 2012). This simplification was appropriate since the fish groups in the base model were not multi-stanza (i.e. subdivided into size or age groups), and size structure of the population was not modelled.
Though the literature often describes respiration as an exponential increase of base metabolic energy demand with temperature, many experimental studies only reported two-to four-fold increases in standard or resting metabolism (Bernreuther et al., 2013;Dalla Via et al., 1998;Johansen and Jones, 2011;Sandersfeld et al., 2017;Schwieterman et al., 2019;.
Studies investigating the relationship between temperature and metabolic costs have included SDA or locomotion, but similar ranges of metabolic increases were reported (Fu et al., 2009;. Therefore, we thought that including resting metabolism, SDA, and activity should only increase energetic losses to metabolism by a factor of eight to ten-fold over biologically relevant temperatures.
The modified Arrhenius equation reported in Blanchard et al. (2012) calculated with their reported parameters was used as the functional form for the thermal response of fish respiration in relationship to temperature (Eq. (4)). The thermal modifier (τ), ranging from zero to ten, is a function of the temperature in degrees kelvin, the Boltzmann constant (k; 8.62x10 -5 eV K -1 ), the activation energy (0.63 eV; similar to values in other studies such as Gillooly et al. 2001, Brown et al. 2004, and a constant (c1=25.55).
The modified Arrhenius equation is the best fit for constraining the metabolic modifier to the range of increases seen in literature for these and similar species.
However, there are no species-specific parameters. Given that any species-specific thermal responses would be clouded in the base model's aggregation to functional groups, we considered the Blanchard et al. (2012) equation to adequately represent a generalized (i.e. non-species specific) fish metabolic thermal response. The respiration thermal response was scaled using similar methods to the consumption response scaling, so that a modifier of τ=1 on the relative respiration curve (RelR; i.e. scaled Blanchard curve) was associated with TB94 (Eq. (5)).
Production rate (P/B) is static, and respiration is solved for in the default programming of Ecopath. Rsim uses a parameter, ActiveRespFrac, to represent the fraction of energy devoted to respiration, which is calculated from the production to consumption ratio and unassimilated food (Aydin et al., 2016). The ActiveRespFrac parameter is carried through the dynamic Rsim simulations. Unassimilated food is assumed to be independent of temperature.
The first step in creating the time series used to force temperature-dependent respiration was to multiply the relative consumption curves by the Ecopath total consumption, equal to Ecopath Q/B multiplied by biomass, to make a curve of total consumption by temperature (Eq. (6)).
The total consumption at TB94 for each of the functional groups was multiplied by the base ActiveRespFrac to get the total respiration at TB94. The relative respiration curve was multiplied by the ratio of the total respiration at TB94 to the relative respiration modifier value at that temperature which produced a total respiration by temperature curve (Eq. (7)).
Dividing the total respiration by total consumption gave ActiveRespFrac by temperature (Eq. (8)).

= ⁄ (8)
In another sensitivity test, we scaled the consumption and respiration curves to the temperatures that informed the original parameters estimated in the Ecopath model to examine the sensitivity of the curves to the implicit thermal parameterization of the  (9)).

Comparison of model versions
Three model versions with two warming scenarios were compared ( Table 1). The

Temperature
There was strong interannual variability in the observed, averaged surface temperatures from 1994-2018; however, there was no apparent warming trend ( Figure   2

Consumption response curves
The thermal response curves for consumption generally reflected the shape as suggested by theory ( Figure 3A). The curves varied between functional groups with piscivorous fish having the steepest increase in consumption below ToptC, while benthivorous fish had the steepest decrease near their thermal maximum. The input species for these curves generally had broad thermal tolerances, and the average temperatures experienced in Narragansett Bay were on the lower end of their tolerance ranges (Supplemental Figure B2.1).
The thermal modifiers of all functional groups increased over 1.0 early in the time series, though within each functional group the modifier only varied by 0.2-0.5 ( Figure 3B). The ending consumption modifier was larger in the high warming scenario for all fish groups, and piscivorous fish had the largest consumption modifier overall (maximum=1.38). The community composition sensitivity test resulted in relative consumption curves with different thermal maxima, but the consumption modifier time series were similar as the curves overlap at the colder average temperatures of the Bay (Supplemental Figure B2.2). The curves created from a single year's observed biomasses were more extreme than the curves created with moderate biomasses of new southern species.

Respiration response curves
The original Blanchard curve had a thermal modifier of 1.0 at 13°C, thus the Blanchard curve and the relative respiration curve scaled to TB94 were similar ( Figure   4A). The three ActiveRespFrac curves by temperature varied by functional group due to the interplay between the respiration and consumption curves ( Figure 4B). The planktivorous fish ActiveRespFrac curve was relatively flat before steeply increasing near maximum temperatures. The benthivorous fish curve showed an increasing trend as waters warm (from 0.54 to 0.74 from 0°C to 20°C). The piscivorous fish ActiveRespFrac decreased to a minimum of energy allocated to metabolism (41%) occurring at 18.4°C before sharply increasing. The temperature-dependent production, shown by the balance of temperature-dependent consumption and respiration, varied in shape by functional group (Supplemental Figures B3.1A-B3.1C). Piscivorous fish had the most defined temperature-production curve, in which the temperature of maximum potential production was slightly cooler than the temperature of maximum consumption. The planktivorous fish production-temperature curve was less pronounced, and there was a very minor curve in the energy available for production for benthivorous fish.
ActiveRespFrac forcing varied between the high and low projected warming scenarios ( Figure 4C) Under the second sensitivity test, the scaling temperature chosen for the relative consumption and respiration curves varied by functional group. The TQB for piscivorous fish was much higher than TB94, which was reflected in the relative consumption curve for that functional group (Supplemental Figure B3.2A). The consumption for piscivorous and benthivorous fish remained below the respective Ecopath starting points for the entirety of the time series (Supplemental Figure   B3.2B). For respiration, the TPB was similar to TB94, so that the ActiveRespFrac by temperature curves differed slightly (Supplemental Figures B3.3A-B3.3B). The ActiveRespFrac time series had similar increasing or decreasing patterns by functional group but were scaled differently relative to the initial ActiveRespFrac (Supplemental Figure B3.3C). In this test, piscivorous fish had respiration costs higher than their baseline in all years.

Comparison of model versions
The inclusion of temperature-dependent fish bioenergetics impacted the modelled biomasses ( Figure 5). For planktivorous and benthivorous fish, the consumption model versions yielded higher ending biomasses than the base model, and the respiration model versions yielded lower biomasses in warming water ( Figure   6). Planktivorous fish and benthivorous fish biomasses were 14.55 g/m 2 and 9.93 g/m 2 in the high warming respiration version compared to 17.03 g/m 2 and 12.46 g/m 2 in the high warming consumption version. Piscivorous fish biomass was highest in the high warming respiration model version (9.25 g/m 2 ) due to their respiration costs decreasing as water temperature increased from its current state. In the model versions, many of the vulnerabilities of strong predation connections for the fish groups were less than 2.0, which corresponds to bottom-up forcing. In the third sensitivity test of changing vulnerabilities, we made these important vulnerabilities more top-down (i.e. closer to 2.0), such that the fish groups were able to further increase their consumption (Supplemental Figures B4.3-B3.5).
Planktivorous fish, in particular, had higher projected biomasses when vulnerabilities were more top-down (Supplemental Figure B4.6). The smaller changes in vulnerabilities between the fitted and 20% change in vulnerabilities yielded smaller changes in biomass (Supplemental Figure B4.7).
The ActiveRespLoss is the parameter adjusted with the new respiration forcing of ForcedActresp. ActiveRespLoss is a function of total consumption, the fraction of energy devoted to respiration, and the forced respiration modifier. All model versions had an ending ActiveRespLoss higher than that of the static Ecopath model ( Table 2).
There were noticeable, but not strong, differences in the 2054 respiration losses, reflecting the variability of the ending respiration modifier and biomasses between model versions.
Ecosim-type models forced into the future can reach an equilibrium state if not forced by varying external drivers Therefore, the total energy inputs and outputs had nearly equilibrated by 2054 in all model versions (Supplemental Table B4.2). Given that fish groups represent just a portion of the biomass of the ecosystem, certain ecosystem metrics, such as catch, were more strongly influenced by changing fish bioenergetics than others. The respiration model versions had the highest total ecosystem respiration and the lowest overall production.

Discussion
We have expanded existing software to make a bioenergetics-based, temperature-dependent food web model and have shown that this new functionality introduced into the flexible Rpath package can impact biomass projections. Climate change elicits complex responses from organisms and ecosystems (Roessig et al., 2004), and our work adds a critical modelling component of temperature-dependent energetic losses. Our methodology can be used in combination with other tools to explore, more completely, the impacts of warming water on marine food webs.

Temperature
The effects of rising temperatures are increasingly being incorporated into

Consumption thermal response
The average thermal response for a functional group can shift as community composition changes, potentially mirroring the changing temperature conditions (Flanagan et al., 2019). Warming waters can restructure marine communities as species move to remain within their preferred thermal habitat (Burrows et al., 2019;Hale et al., 2017;Nicolas et al., 2011). Our sensitivity test on changing community composition showed that large interannual variability in relative species abundance can have as strong an impact on the consumption thermal response as the introduction of new species with warmer tolerances if those new species begin at relatively small biomasses. However, we would expect the warm-water species to become more dominant over time as the environment becomes more favorable for them. Such a shift in species has been seen in Narragansett Bay over the last few decades (Collie et al., 2008;Oviatt et al., 2003). Changes in community structure could be explored further in future models that split warm and cold-water species into distinct functional groups.

Respiration thermal response
The Blanchard et al. (2012) parameterization of the Arrhenius equation yielded a reasonable description of changing metabolic costs for the fishes in our ecosystem.
Ideally, we would have had species-specific responses, but since there were few respiration data available for the species in our ecosystem, we would have used assumed values regardless. Some studies report a temperature that corresponds to a maximum resting metabolism beyond which metabolic demand decreases (Bernreuther et al., 2013;MacIsaac et al., 1997;Schulte, 2015), but this decrease is not always seen (Giacomin et al., 2017;McKenzie et al., 2016; Dabrowski, 1985;Priede, 1985;Sun et al., 2006). We strongly recommend ground truthing ActiveRespFrac values in the balancing step when building future EwE and Rpath models to be used to examine thermal drivers and bioenergetic questions. Previous work has also found that temperature altered growth rates as a result of the balance between energy inputs and outputs (Cotton et al., 2003;Gaylord et al., 2003;Present and Conover, 1992). The production curves varied between the fish groups, and the aggregation of multiple fish species with differing data quality into a single functional group response was likely responsible for any deviations in curve shape from theory.
Ecopath models are implicitly parameterized for ambient temperatures, and when the temperature changes, so will the vital parameters. The initial piscivorous fish life history parameters (i.e. Q/B and natural mortality) from Fishbase were generally estimated from warmer temperature environments which can result in higher consumption and metabolism (Froese and Pauly, 2019). Creating Ecopath models can necessitate borrowing parameter values from similar species and geographic regions, but users should consider inconsistencies in environmental conditions if they hope to include temperature impacts in their modeled ecosystem.

Comparison of model versions
Biomass projections between the model versions were noticeably different, though these differences were small due to the consistent forcing of phytoplankton and fishing mortality stabilizing the outputs. The fish biomass projections were similar for 1994-2018 and diverged beyond 2018 as the temperature increased. The consumption versions had the highest biomasses for benthivorous and planktivorous fish reflecting the increased energy intake and assuming that all of the energy was available for production. The respiration model versions had the lowest benthivorous and planktivorous fish biomasses, a result consistent with increased energetic demands of warmer waters .
Ignoring the shifting bioenergetic balance when modelling fish biomasses in response to climate change may yield optimistic forecasts for fisheries. Changing productivity of fish populations can potentially lead to lower sustainable fisheries yields (Free et al., 2019), and different management approaches may need to be considered when examining a warming ecosystem (Serpetti et al., 2017). As top predators, piscivorous fish likely had the smallest biomass differences between the respiration and consumption model versions because they were both externally forced with changing energetic demands, and their prey (i.e. the other fish groups) were also affected. In the respiration model versions, the respiration costs for piscivorous fish were decreasing, and the group might have had higher biomasses if their prey fish had not declined due to their own increasing metabolic costs. Trophic amplification of environmental impacts has been shown in other models of warming systems, where predators showed greater declines than forced declines in primary producers because the intermediate groups were also affected (Kwiatkowski et al., 2019;Lotze et al., 2019). Even limiting variable forcing to the fish groups had trickle-down effects on other groups of the ecosystem. While the changes in the biomass of other functional groups were more minor than the fish groups, they were still noticeable which highlights the importance for bioenergetic considerations when projecting food web impacts to climate change.
The vulnerabilities estimated for our fish groups were generally bottom-up, so the predation mortality that these fish groups can exert on their prey was limited even as fish biomass increases (Christensen and Walters, 2004). Our third sensitivity tests illustrated the influence of vulnerability parameters on Ecosim and Rsim biomass projections, which has been described by others SEDAR, 2020). Because higher (i.e more top-down) vulnerabilities led to higher consumption and biomasses, top-down ecosystems may be more strongly affected by changing predator bioenergetics.
We chose to not refit the model versions (i.e. re-estimate vulnerabilities) after including the temperature drivers to isolate the impacts of changing bioenergetics.
Therefore, while biomass projections differed, we did not necessarily achieve a better model fit with the inclusion of temperature-dependent fish bioenergetics. A model with vulnerabilities estimated by fitting procedures after adding the thermal responses may better match observed biomass as has been seen in other EwE models (Bentley et al., 2020(Bentley et al., , 2017. However, our focus was on developing new metabolic functionality in Rpath and understanding the scale at which temperature-dependent fish bioenergetics can impact predicted biomasses rather than using bioenergetic theory to explain observed trends.

Considerations for future models
This model required a vast amount of data, but data for basic bioenergetic parameters do not always exist, particularly as studies may be biased toward species of economic importance or those able to thrive in laboratory conditions . Physiological studies seem to emphasize the nuance of fish response to temperature. While advanced bioenergetic models may be able to capture this nuance, We modeled simple bioenergetic responses to temperature based on established bioenergetic principles. In reality, the factors influencing consumption and metabolism are much more complex. Thermal responses can differ between individuals as well as populations, and responses can be altered by the presence of simultaneous stressors (Farrell, 2016;Kroeker et al., 2013;Present and Conover, 1992). A more detailed base model may address the stage-specific thermal tolerances and the role of body size in metabolic performance (Hare et al., 2010;Sibly et al., 2012), but this would require additional size-or age-specific bioenergetic data for each species. We did not include acclimation effects, though the previous thermal exposure of organisms can affect their response to temperature, nor the rate of temperature change, which can alter an organism's response (Morgan et al., 2018;Pinsky et al., 2020). Our work also excludes the feedbacks between temperature, reproduction, and recruitment processes (Conover and Kynard, 1981;James, 2020;Johnston et al., 1998;Pankhurst and Munday, 2011), which may be included in future models with increased use of multi-stanza groups. Finally, individual organism behavioral responses such as behavioral thermoregulation or changes in risk taking behavior were not included (Nagelkerken and Munday, 2016; Neubauer and Andersen, 2019).
Our work could be expanded in the future to include thermal responses experienced by other non-fish functional groups. Climate change has been shown to influence community composition, abundance, and timing of plankton blooms (Lawrence and Menden-Deuer, 2012;Smith et al., 2010;Sullivan et al., 2001), which could have significant feedbacks on Narragansett Bay energy flow (Monaco and Ulanowicz, 1997). Reduced primary production due to climate change could result in reduced fisheries catch if lower level production can no longer support high predator abundance (Brown et al., 2010;Cheung et al., 2011;Johansen et al., 2015). Fishing and other human activities are also considered significant drivers of ecosystem indicators (Link et al., 2010). More confidence and precision in future biomass estimates could be achieved through the inclusion of the thermal responses of the lower trophic level groups and additional varying top-down and bottom-up dynamics.

Conclusion
EwE has most often been used to address fisheries harvest questions (Pauly et al., 2000); our work builds on previous applications of environmentally driven food web models used to explore climate change impacts. In our Narragansett Bay case study, we have demonstrated that the inclusion of temperature-dependent bioenergetic drivers into existing Rpath and Rsim models can alter projections of biomass and energy flow. Ecosystem health and sustainability goals can be better achieved by integrating principles of conservation physiology into multispecies models to gain improved resolution on how ecosystems will respond to environmental pressures . Productivity of Narragansett Bay, or similar ecosystems, may be compromised in a warmer future if fish or other ectotherms devote greater amounts of energy towards meeting metabolic demands. Our novel methodology is intended to serve as an example to address such questions in other warming marine ecosystems (Pershing et al., 2015) or lake environments (Adrian et al., 2009) in which there may limited ability for fish to seek alternate temperatures. The new functionality for temperature-dependent respiration in the Rpath package can be combined with other thermally-driven model components to provide the most comprehensive predictions of how a food web will respond to climate change.

Acknowledgments
This work was supported by Rhode Island National Science Foundation Established Program to Stimulate Competitive Research Grant #OIA-165522. We acknowledge the World Climate Research Programme, which, through its Working Group on Coupled Modelling, coordinated and promoted CMIP6. We thank the climate modeling groups for producing and making available their model output, the Earth System Grid Federation (ESGF) for archiving the data and providing access, and the multiple funding agencies who support CMIP6 and ESGF. We thank those that assisted with the development of this model and for their edits on this manuscript including J. Langan, C.
Suckling, and R. Bell. This work is a contribution of the Rhode Island Marine Fisheries Institute.

References
Adrian           (Froese and Pauly, 2019) length at first maturity was used to make a determination in absence of other data.

Maximum temperature
The maximum temperature was rounded to a whole degree using the ceiling function to ensure the rounded maximum temperature encompassed the reported temperature.
Three main data sources were considered, and the highest value for each species were chosen as the final maximum temperature. The three data sources were  ). 6. In the absence of other data, the temperature of maximum consumption is assumed to be 90% of the maximum temperature. This is based on the typical shape of the Kitchell curve as it is described as a skewed curve with a steep drop off. That shape can only be achieved with a temperature of maximum consumption near the maximum temperature.
Supplemental The consumption curve created with that high Q10 did not match the shape as specified in bioenergetic theory. The Q10 value was lowered by choosing only the estimates of feeding from 6°C to 10°C.
Species-specific biomass time series were used to create a weighted average of the thermal response curves by functional groups (Supplemental Table A3.4). The time series were the same as those used in Innes-Gold et al (2020).
Supplemental Table A3. 4: 1994-1998 averages of species-specific biomasses were used to weight the thermal parameters. Data were taken from the GSO and DEM trawls, with more detail found in . Biomasses in g/m 2 have been rounded to four decimal points.  Table A4.1: Southern species used to assess sensitivity of the consumption curve to changing community composition. Species were chosen as those that have been caught in Narragansett Bay in the last 15 years according to the GSO trawl data and are currently found in the Chesapeake Bay Jung and Houde, 2003). Functional group was assigned after consulting the diet as reported in Bowman et al. (2000). All have warmer thermal maxima than the species of the original model. Maximum temperature was taken as the extreme value of the temperature of catch from the GSO fish trawl and the website Aquamaps (accessed March 2021). The temperature of optimum consumption was assumed to be 90% of the maximum temperature, and Q10 was assumed to be 2.3. The biomass values were chosen as the median biomass of the original species that made up the functional group, so the species was represented similarly the others of the functional group.  (3), except the Kitchell modifier was evaluated at the TQB of each functional group instead of TB94. The respiration thermal response was scaled similarly to what was done for the consumption response, so that a modifier of 1.0 on the scaled respiration curve (i.e. Blanchard curve) was associated with the temperature that informed the Ecopath P/B parameter (TPB). In EwE, the P/B parameter is estimated as total mortality (Z) which is the sum of natural mortality (M) and fishing mortality (F) . The natural mortality input literature generally reported temperatures, and the temperature of fishing mortality was assumed to be the 1994-1998 yearly averaged temperature of Narragansett Bay as the F values were calculating from Bay catches. In cases for which temperature of natural mortality was not reported, it was assumed to be the average of the 5 th and 95 th percentiles for temperature reported in Aquamaps . The temperature of natural mortality (TM) and fishing mortalities (TF) were averaged to determine TPB. The relative respiration was calculated according to Eq. (5) except that the original Blanchard modifier was evaluated at the mean of TM and TF instead of TB94.

Species
The temperature of Ecopath (TEco) for each functional group was considered to be the average TQB and TPB, described earlier. Total respiration was set so that, at TEco, total consumption divided by total respiration was equal to the original ActiveRespFrac of the base model version. Total respiration was calculated according to Eq. (7) except that the total consumption and relative biomass curves were evaluated at TEco instead of TB94. All temperatures by functional group are listed in Supplemental