Big in the benthos: Future change of seafloor community biomass in a global, body size‐resolved model

Abstract Deep‐water benthic communities in the ocean are almost wholly dependent on near‐surface pelagic ecosystems for their supply of energy and material resources. Primary production in sunlit surface waters is channelled through complex food webs that extensively recycle organic material, but lose a fraction as particulate organic carbon (POC) that sinks into the ocean interior. This exported production is further rarefied by microbial breakdown in the abyssal ocean, but a residual ultimately drives diverse assemblages of seafloor heterotrophs. Advances have led to an understanding of the importance of size (body mass) in structuring these communities. Here we force a size‐resolved benthic biomass model, BORIS, using seafloor POC flux from a coupled ocean‐biogeochemistry model, NEMO‐MEDUSA, to investigate global patterns in benthic biomass. BORIS resolves 16 size classes of metazoans, successively doubling in mass from approximately 1 μg to 28 mg. Simulations find a wide range of seasonal responses to differing patterns of POC forcing, with both a decline in seasonal variability, and an increase in peak lag times with increasing body size. However, the dominant factor for modelled benthic communities is the integrated magnitude of POC reaching the seafloor rather than its seasonal pattern. Scenarios of POC forcing under climate change and ocean acidification are then applied to investigate how benthic communities may change under different future conditions. Against a backdrop of falling surface primary production (−6.1%), and driven by changes in pelagic remineralization with depth, results show that while benthic communities in shallow seas generally show higher biomass in a warmed world (+3.2%), deep‐sea communities experience a substantial decline (−32%) under a high greenhouse gas emissions scenario. Our results underscore the importance for benthic ecology of reducing uncertainty in the magnitude and seasonality of seafloor POC fluxes, as well as the importance of studying a broader range of seafloor environments for future model development.


| INTRODUCTION
Of the particulate organic carbon (POC) fixed annually in the surface ocean by primary producers (% 45 Pg C year À1 ; Behrenfeld & Falkowski, 1997), almost all is remineralized back to dissolved inorganic carbon (DIC) within a short time period (months to years). However, mediated by a complex community of benthic organisms, a small fraction of this POC (0.05 Pg C year À1 ; Hain, Sigman, & Haug, 2014) is sequestered by seafloor burial, and represents a significant flux on geological timescales (Mawbey & Lear, 2013). This same community also remineralizes POC, turning it back to DIC and nutrients that eventually resupply productive surface communities (Dunne, Sarmiento, & Gnanadesikan, 2007). Consequently, understanding the role that benthic communities play in the biogeochemical cycles of the ocean, and estimating how this role may change into the future, is of considerable interest (e.g. Kriest & Oschlies, 2013;Moodley et al., 2011).
Benthic communities are an important source of commercially exploited seafood, both for direct human consumption, and for other harvested species. They are also often highly biodiverse, and in the deep-sea, may be relatively pristine compared to other marine habitats more directly impacted by human activities (Jackson et al., 2001;Lotze et al., 2006;Ramirez-Llodra et al., 2011).
Advances in deep-sea photography (e.g. Bett, 2003), sampling (e.g. Gooday, Bett, Shires, & Lambshead, 1998) and even in situ experimentation (e.g. Jeffreys et al., 2013;Main et al., 2015;Nomaki, Heinz, Nakatsuka, Shimanaga, & Kitazato, 2005;Witte, Aberle, Sand, & Wenzh€ ofer, 2003) have significantly increased understanding of how these communities are organized (Danovaro, Snelgrove, & Tyler, 2014). Evidence suggests that food web complexity and functional groupings of organisms may be of secondary importance to body size in controlling energy flow through these communities (Blanchard, Heneghan, Everett, Trebilco, & Richardson, 2017;Brown, Gillooly, Allen, Savage, & West, 2004;Dickie, Kerr, & Boudreau, 1987;Peters, 1983). As body mass is readily quantified and well correlated with metabolic processes (e.g. Brey, 2010) as well as community biomass and abundance (e.g. Hildrew, Raffaelli, & Edmonds-Brown, 2007), this is attractive on both practical and theoretical grounds. However, there is still considerable debate among biologists and ecologists as to the causes, and universality, of the mass-scaling of metabolism (e.g. Glazier, 2005Glazier, , 2010Hildrew et al., 2007;Hunt & Roy, 2006). Nevertheless, the continuing development of the Metabolic Theory of Ecology (McClain, Allen, Tittensor, & Rex, 2012;Schramski, Dell, Grady, Sibly, & Brown, 2015) demonstrates the practical value of size-based assessments of ecosystem function. Allometry has proven a useful concept in both terrestrial and aquatic ecology (e.g. Hildrew et al., 2007;Schmidt-Nielsen, 1984), and its application to benthic communities, and marine systems in general, has clear potential (e.g. Blanchard et al., 2009Blanchard et al., , 2017. As a preliminary step towards this, Kelly-Gerreyn et al. (2014) introduced an allometry-based model of the seafloor metazoan community that reproduced the distribution of biomass and abundance at three contrasting locations: (i) Faroe-Shetland Channel, North Atlantic deep-sea; (ii) Fladen Ground, North Sea continental shelf; and (iii) Oman Margin, Arabian Sea continental slope (see section 2.1.1 for more details). They employed high-quality size-resolved observational data on the meio-and macrobenthos from these three sites. Their model was tuned for these locations individually and subsequently retuned to all three sites simultaneously to provide a more broadly applicable, or "unified'' parameterization (Ichino et al., 2015; see section 2.1.1 and Appendix S1). The three study locations are highly contrasting in terms of water depth, habitat temperature, and the magnitude and seasonality of productivity. Consequently, Kelly-Gerreyn et al. (2014) consequently suggested that the model was likely to be robust and of broad application to marine benthic communities. (IPCC) found that, under a high greenhouse gas emissions scenario (RCP 8.5; see section 2.2.1), primary production changed by an average of À8.1%, with falls in export production being even greater, ranging from À7% to À16% (Bopp et al., 2013). Similar results were found by Yool, Popova, Coward, Bernie, and Anderson (2013), who inferred a link with ocean acidification (OA), which led to drops in primary production (À6.3%) translating to much larger declines in export flux to the deep ocean (À40.7% at 1000 m water depth).
Using simulated changes in export flux and an empirical seafloor biomass model, Jones et al. (2014) estimated that an average change in seafloor POC flux of À11.4% translated to a decline in total benthic biomass of À5.2% at the global scale. These changes suggest that benthic communities will experience substantial impact in the near future.
In addition to changes in the food supply available to them, although not examined here, it is anticipated that benthic communities will experience other stresses. Temperature has been already been identified as a key stressor from palaeological records (Yasuhara, Okahashi, Cronin, Rasmussen, & Hunt, 2014), especially where it is naturally low or high and organisms are most vulnerable to change (Yasuhara & Danovaro, 2016). More broadly, much as with pelagic communities (Gruber, 2011;Popova et al., 2016), benthic communities are also vulnerable to additional anthropogenically driven stressors, including acidification, deoxygenation and contamination (Levin & Le Bris, 2015). Exposure to stressors such as temperature and acidification would likely be strongly depth-dependent because of their surface sources, but the other stressors may YOOL ET AL.
| 3555 be communicated more directly to the deep ocean through the biological pump (deoxygenation) or gravitational sinking (contamination).
Here we take Kelly-Gerreyn et al. (2014)'s allometry-based model of the benthic community and force it with POC fluxes derived from the pelagic ecosystem model of Yool, Popova, Coward, et al. (2013). The model's unperturbed behaviour is examined, both in terms of global geographical variability and the role of temporal variability (especially seasonality) in POC flux. How these features may change into the future is then investigated using two end-member IPCC AR5 scenarios for the 21st century (RCPs 2.6 and 8.5). We focus on the outcomes for integrated seafloor biomass and its distribution across the modelled size classes, and how both of these properties vary in space and time.

| BORIS model
The benthic community is modelled using the Benthic Organisms Resolved In Size v1.0 (BORIS-1; henceforth BORIS) model (Kelly-Gerreyn et al., 2014), 1 which represents the biomasses of seafloor metazoans in the meio-to macrofaunal size range (0.9 lg wet wt to 30 mg wet wt; see Table S1). This size range was selected based on the availability of high-quality size-resolved field data. POC flux consumption and subsequent respiration by smaller and larger organisms are represented implicitly via the f other parameter (see Appendix S1). As illustrated in Figure 1a, the model is driven by the POC flux reaching the seafloor which, after a fixed fraction is consumed by the implicit respiration of other members of the benthos, enters a detrital reservoir that is accessed by the modelled metazoans, both as a source of food (ingestion) and a sink for losses (mortality, defecation). Ultimately, all of the POC that arrives at the seafloor is consumed and respired; that is, burial sequestration is assumed to be minimal, and the model provides a biomass size distribution of the modelled metazoan community. Long-term observations of POC flux to the seabed and corresponding sediment community respiration in the deep ocean suggest that this is a reasonable simplification (Smith & Kaufmann, 1999;Smith et al., 2009). Similarly, global-scale assessments of the deep ocean (>2,000 m water depth) POC burial flux suggest a value of only 3% of seafloor POC flux (Dunne et al., 2007).
Appendix S1 presents a detailed description of the structure, parameterization and evaluation of BORIS. The computational implementation of BORIS is described in Appendix S1.

| Default performance
In Kelly-Gerreyn et al. (2014), the eight model parameters in BORIS were tuned using a microgenetic algorithm approach (Ward, Friedrichs, Anderson, & Oschlies, 2010). This was carried out separately for each of the three field locations. The procedure is fully described in Kelly-Gerreyn et al. (2014) and summarized in Appendix S1.
The field data used to calibrate BORIS came from: Faroe-Shet-  Table S2 lists the values of the parameters derived by Kelly-Gerreyn et al. (2014) for the three sites, together with the ranges from which they were drawn. Also listed is a "unified'' parameter set (column All) that was produced using the same approach as in Kelly-Gerreyn et al. (2014), but fitting the model to all three sites simultaneously.
This parameter set was previously employed by Ichino et al. (2015), and is used throughout this study. In Kelly-Gerreyn et al. (2014), the model was not explicitly given a name, but we adopt the moniker BORIS here to facilitate reference and discussion. Yool, Popova, Coward, et al., 2013). POC flux to the seafloor comprises two separate size classes of detrital particles, the balance of which means that shallow regions are typically dominated by "small,'' slow-sinking particles, while deep regions are dominated by "large,''
Appendix S2 presents more details concerning the structure, parameterization and validation of NEMO-MEDUSA.

| Forcing scenarios
We In terms of POC fluxes to the seafloor, Yool, Popova, Coward, et al. (2013) found that these declined into the future in parallel with the decline in primary production driven by increasing ocean stratification and decreasing nutrient availability. However, while the decline in global primary production is modest (À6%), the corresponding decline in POC flux to the deep ocean is greater (À40%; 1,000 m water depth), as a result of the role of ocean acidification. Driven by oceanic uptake of anthropogenic CO 2 , this both shoals the modelled calcite compensation depth and decreases NEMO-MEDUSA's production of biogenic calcium carbonate, with the latter the dominant factor for POC fluxes.
Appendix S2 presents more details concerning these scenarios and the forcing output.
These simulations were started from analytically derived initial conditions (see Appendix S1), and then spun-up for a period of This analysis was then extended using historical  and future   In all simulations of BORIS, the seafloor POC fluxes provided by NEMO-MEDUSA were first amended by subtracting a fixed and geographically uniform proportion of POC-the parameter f other . As noted above, this is designed to account for the fraction consumed by unmodelled seafloor organisms (e.g. microbes and megabenthos).
The remaining POC entered the seafloor detritus pool at which point it became available for BORIS's modelled metazoan size classes. This parameterization is discussed further later, and details concerning it can be found in Appendix S1. As noted earlier, organisms outside the size range considered by BORIS (larger and smaller) consume a portion of the POC flux, and the feeding relationship of all of these organisms on seafloor POC is much more complex than the simple transfer to respiration represented here.

| RESULTS
In the following, BORIS only implicitly represents both smaller and larger size classes that contribute to complete seafloor biomass.
Consequently, the term "total biomass'' is used in the presentation of results from BORIS to refer specifically to the total biomass of the modelled size classes resolved by BORIS and not to the complete total seafloor biomass of living organisms. YOOL ET AL.
| 3557 Table 1 presents the observed and modelled seafloor POC fluxes for the three sites alongside both those observed, and those from the optimized BORIS model (per Table S2). The input POC flux that is available to the modelled size classes of BORIS, Q, is a fraction of the total flux reaching the seafloor, and the observed fluxes in

| Global patterns
Although the three calibration sites examined above vary in both total POC flux and in seasonal pattern, they represent only a fraction of the global range in these properties. Figure    The general pattern in average temporal lags between peak annual seafloor POC flux and seafloor detritus (Fig. S5b) or modelled total biomass (Fig. S5c) is for the lag to increase with decreasing POC flux, with a tendency for shorter lags where CoV is higher (more seasonal). Fig. S6 presents the corresponding lags for size classes 1 and 16, and illustrates a broadly consistent pattern of greater lags in the larger size classes (see also Fig. S7).

| Future change
To investigate the impact of changes in seafloor POC flux for benthic communities, Figures 4, S8 and Table 2 detail the change in globally integrated seafloor biomass across the 21st century for the RCP 2.6 and RCP 8.5 scenarios described earlier (Rogelj et al., 2012). Table 2 reports globally integrated primary production, its translation to export production and the flux of POC to the seafloor, and the resulting modelled total biomass at a series of water depth intervals from shallow seas (to 0.1 km) to the full depth of the ocean (to F I G U R E 2 NEMO-MEDUSA simulated seafloor POC flux, as (a) mean and (b) coefficient of variation (CoV; based on monthly means to illustrate seasonality) F I G U R E 3 (a) Mean annual field of total modelled seafloor biomass (Fig. S3 shows the corresponding field of seafloor detritus).
(b) The relationship between mean annual POC and mean annual seafloor detritus (R; black) and biomass (classes 1-16; colours). Note that POC flux is expressed as g wet wt m À2 day À1 , and that logarithmic scales are used throughout scenario experiences a 6.1% decline by the 2090s, control of emissions in RCP 2.6 results in a small increase in global productivity by the end of the century (in part related to slightly elevated temperatures enhancing phytoplankton growth rates). Export production shows a small decline of À1.3% under RCP 2.6 (in part related to slightly elevated temperatures enhancing remineralization rates), but a much larger decline of À11.4% under RCP 8.5's more extreme change. Integrated to the global scale, POC fluxes to the seafloor are increased for RCP 2.6 (+8.5%), while under RCP 8.5 there is a decline (À3.9%). This seeming disparity with export production stems from the seafloor POC flux being biased towards shallow regions (where a greater fraction of export survives to the seafloor), and the exclusion of regions with depths <100 m from the reported export production.
To separate the impact of changes in seafloor POC flux and water depth, Table 2 divides the seafloor into seven water depth bands. For the shallow seas (<100 m water depth), both scenarios result in an increase in seafloor biomass, largely driven by warmer conditions that increase growth rates in shallow regions (Yool, Popova, Coward, et al. 2013). For both scenarios, this increase declines and then reverses (by 200 m water depth under RCP 8.5) with increasing water depth. By 5,000 m water depth, seafloor biomass is lower for RCP 2.6 (À7.0%) and much lower for RCP 8.5 (À32.0%).
Integrating globally, the total modelled seafloor biomass in BORIS declines slightly under RCP 2.6 (À1.1%), contrary to the corresponding increases in both primary production and seafloor POC flux.
To distinguish the influences of changing primary productivity and water depth, Table 3  or greater) and modelled total biomass (À32.0% or greater). The relationship between change in seafloor POC flux and change in modelled total biomass is strongly linear with a slope of approximately 0.6, that is biomass changes less than the POC flux. This contrasts with the pelagic situation of NEMO-MEDUSA where a decline of primary production of À6.1% was much more closely paralleled by a corresponding decline of À5.7% in the biomass of surface phytoplankton (Yool, Popova, Coward, et al. 2013). Table 2 and Fig. S9 indicate how changes in primary production, export production, seafloor POC flux and seafloor biomass are correlated. In the case of export and primary production (Fig. S9a), this is almost linear, with a slight bias towards decreased export; that is, a À10% decline in primary production leads to a large decline in export, while a 10% increase in primary production leads to a slightly smaller increase in export. Similar, but less well-correlated and much less linear, relationships are apparent between seafloor POC flux and both export (   Temporal variability can also be considered in terms of the approach to equilibrium of body mass distribution. Through idealized perturbations of the seafloor detritus pool (half and double), Fig. S11 shows that the time to equilibrium for BORIS increases with decreasing seafloor POC flux, that is communities with lower total biomass take longer to return to equilibrium, and that it may take a  ger disequilibrium although the biomass in these regions is much lower (Fig. S13a). At the end of the century, there is a greater positive disequilibrium in most areas, and noticeably larger areas of positive disequilibrium in the gyres (Fig. S13b). This general shift to more positive disequilibrium stems from a combination of declining ocean productivity and the system lags shown in Fig. S11, especially for the low biomass/slow equilibration gyre regions. Fig. S13c, d show the most "out of equilibrium'' state variable of BORIS for the 1990s and the 2090s.
Broadly, the global benthos is split between locations with excess seafloor detritus, that is more than at equilibrium, and those where the largest size class is most out of equilibrium. By the 2090s, the balance of this situation has shifted somewhat, such that the largest size class is most out of equilibrium over a larger portion of the world.
In general, intra-and interannual variations in seafloor POC flux introduce deviations, but as the three calibration sites illustrate, the dynamical simulation time-averages the equilibrium state quite closely (Fig. S12b). Consequently, the future behaviour of BORIS can be accurately estimated using time-averaged seafloor POC fluxes and steady-state calculations, at least for areas with non-negligible POC fluxes and shorter equilibration times.
A key result from this study is a forecast decline in benthic biomass substantially greater than the corresponding decline in surface ocean productivity that drives it (see also Yool, Popova, Coward, et al. 2013). Jones et al. (2014) also examined the fate of benthic biomass in a multimodel study, finding a decline in productivity and near-seafloor POC flux (2000s ! 2090s) of more than À11% that resulted in a decline in T A B L E 3 Seafloor area (10 6 km 2 ), POC flux to the seafloor (mg C m À2 day À1 ) and modelled total biomass (g wet wt m À2 ) listed at the global scale and for seven water depth bands for the 1990s and 2090s under scenarios RCP 2.6 and RCP 8.5. For both scenarios, the percentage change between the 1990s and 2090s is indicated in brackets seafloor biomass of more than À5%. Our results appear more marked, so to make a direct comparison, we have forced the empirical model of Jones et al. (2014) with the same NEMO-MEDUSA output used in the present study (Appendix S4). Table 4 indicates that when subject to common forcing regimes, the two models produce very similar outcomes. Comparing the broadly similar body size ranges of total biomass in BORIS with macrofauna biomass in Jones et al. (2014), the respective changes are À6.2% and À7.3% under RCP 2.6 and À31% and À37% under RCP 8.5. These closely matched results derived from substantially different approaches point to the details of the POC flux forcing as being more important than the type of the benthic model being forced by that flux.
This study couples two models, each with their own limitations.
Use of additional models of POC forcing could increase confidence in results; however, the performance of NEMO-MEDUSA has been extensively validated Yool, Popova, Coward, et al., 2013), and the present results suggest that the behaviour of BORIS is a relatively simple function of seafloor POC flux. Only monthly average POC flux forcing has been applied in the present case, and it is possible that shorter-term variability may operate in the field (Smith et al., 2014;Witte et al., 2003).
The BORIS model makes a number of assumptions and simplifications that influence its representation of the benthos and their ecology. As currently configured, BORIS omits organisms outside the reliable range of the calibration data. BORIS effectively represents seafloor microbes as "external'' to the modelled system; that is, they are ignored after they take a portion of the POC flux, and this overlooks the potential role of meio-and macrobenthos in preconditioning detrital substrates (Rowe & Deming, 2011). A practical consequence is the requirement for the f other parameter, which makes specific assumptions and must be separately derived. It is possible to extrapolate the optimized model to both larger and smaller size classes (Ichino et al., 2015), but this was not carried out here in order to maintain model traceability with Kelly-Gerreyn et al. (2014).
For example, the Metabolic Theory of Ecology (Brown et al., 2004) would predict an inverse relationship between habitat temperature and standing stock biomass, with a potential effect factor of 2 between extremes of present-day deep-sea temperatures. Similarly, recent work has highlighted the role of increasing temperature in reductions of body size (Reuman, Holt, & Yvon-Durocher, 2014 (Dunlop et al., 2016;FitzGeorge-Balfour, Billett, Wolff, Thompson, & Tyler, 2010). However, the flux of material from the surface ocean will also include organic matter that has been extensively reworked in the surface and/or mesopelagic realms, which will be of much lower quality (Valls et al., 2014 | 3563 the first-order seafloor POC flux. As such, it overlooks nutritional factors that may play a role in real systems; however, knowledge of the significance of such aspects for benthic ecology is incomplete. Despite these limitations, BORIS represents a dynamic, timedependent alternative to empirical models such as Wei et al. (2010), and several avenues for future refinement are clear: extending the size range of organisms; including explicitly resolving microbes and their role in remineralization; ecological relationships between metazoans and microbes (microbial gardening; Mayor, Sanders, Giering, & Anderson, 2014) and between metazoans (predation).
A key driver of future work will be the assimilation of additional observational data sets, ideally from diverse sites that differ from those used in the initial development of BORIS. Our results on equilibration times would suggest, in particular, examination of sites with very low seafloor POC fluxes to provide end members that can constrain model behaviour under conditions of extreme oligotrophy.
Additionally, similar time-lagged responses to climatic and upper ocean processes have been detected or inferred in abyssal benthic systems (Laguionie-Marchais, Paterson, Bett, Smith, & Ruhl, 2016;Ruhl & Smith, 2004). Incorporation of data from such sites will increase confidence when extrapolating to future situations in which local conditions have significantly departed from current. Such sites may also highlight environmental factors that are currently omitted from BORIS and which may drive its further development. Similarly, the assimilation of new seafloor POC flux data sets also represents a critical future avenue, especially where high temporal frequency sampling is achieved (Smith et al., 2014).
From a model-focused perspective, a potential extension to the work reported here would be to utilize BORIS in simulations forced with the output from a broader range of pelagic ecosystem models.
For instance, from either the existing CMIP5 database (Jones et al., 2014), or those models participating in the upcoming CMIP6 exercise. In addition to evaluating BORIS under a broader range of potential present-day and future conditions, exposure to different models may provide an impetus to consider factors beyond simple bulk POC flux, that some of these models, including NEMO-MEDUSA, can provide.
In summary, here we apply a body mass allometry-based model both for which benthic components should be monitored, and for the frequency at which that monitoring should occur. Under future climate change scenarios, modelled benthic biomass was found to decrease to a greater degree than that of surface ocean productivity (À18% as compared to À6%), with deep-water communities experiencing greater declines than those in shallow seas. This is in accordance with patterns of seafloor POC fluxes in NEMO-MEDUSA, which are increasingly attenuated with water depth by future change, and with the dominant role of seafloor POC flux magnitude in BORIS. We note that absent forcing factors (temperature, oxygen, pH), and our simplified ecological assumptions (implicit microbes, trophic relations), represent important aspects for future development of BORIS, with the acquisition of observations from a broader range of seafloor environments of key importance.

ACKNOWLEDG EMENTS
The authors are grateful for the input and support provided by