A Network‐Scale Modeling Framework for Stream Metabolism, Ecosystem Efficiency, and Their Response to Climate Change

Climate change and the predicted warmer temperatures and more extreme hydrological regimes could affect freshwater ecosystems and their energy pathways. To appreciate the complex spatial and temporal interactions of carbon cycling in flowing waters, ecosystem metabolism (gross primary production [GPP] and ecosystem respiration [ER]) must be resolved at the scale of an entire river network. Here, we propose a meta‐ecosystem framework that couples light and temperature regimes with a reach‐scale ecosystem model and integrates network structure, catchment land cover, and the hydrologic regime. The model simulates the distributed functioning of dissolved and particulate organic carbon, autotrophic biomass, and thus ecosystem metabolism, and reproduces fairly well the metabolic regimes observed in 12 reaches of the Ybbs River network, Austria. Results show that the annual network–scale metabolism was heterotrophic, yet with a clear peak of autotrophy in spring. Autochthonous energy sources contributed 43% of the total ER. We further investigated the effect of altered thermal and hydrologic regimes on metabolism and ecosystem efficiency. We predicted that an increase of 2.5°C in average stream water temperature could boost ER and GPP by 31% (24%–57%) and 28% (5%–57%), respectively. The effect of flashier hydrologic regimes is more complex and depends on autotrophic biomass density. The analysis shows the complex interactions between environmental conditions and biota in shaping stream metabolism and highlights the existing knowledge gaps for reliable predictions of the effects of climate change in these ecosystems.

One of the biggest challenges faced by hydrologists and stream ecologists lies in the attempt to integrate and resolve ecosystem dynamics at the scale of entire river networks. The advent of affordable sensors has boosted the availability of reach-scale time series of GPP and ER estimated from dissolved oxygen (DO) time series (Appling et al., 2018;Beaulieu et al., 2013;Diamond et al., 2021;Hall & Beaulieu, 2013), greatly advancing our understanding of drivers of local stream metabolism (e.g., Bernhardt et al., 2018;Segatto et al., 2020;Ulseth et al., 2018). However, few studies have so far attempted to upscale GPP and ER from stream reaches to an entire river network (see Rodríguez-Castillo et al. [2019] for the spatial distribution of metabolic metrics and Segatto et al. [2021] for network-scale ecosystem regimes both in space and time). We argue that the lack of suitable modeling approaches and distributed data has hampered the development of river metabolism models integrating both river network structure and ecosystem functioning at the scale of entire stream networks. Furthermore, scientific advancement in this direction is required in order to draw conclusions on catchment, regional, and global scales carbon fluxes, to unveil emerging properties at the scale of entire stream networks, and to investigate the effect of climate change.
In this study, we propose a spatially distributed, process-based modeling framework to study and predict metabolic regimes over an entire stream network. The river network is viewed as a meta-ecosystem where a set of ecosystems are spatially connected by flows of energy and materials (Loreau et al., 2003). The model takes advantage of recent results detailed in Segatto et al. (2020Segatto et al. ( , 2021 and couples the thermal and light regimes reconstructed via machine learning with a reach-scale ecosystem model. The model integrates information on network structure, catchment vegetation, and hydrological regime to simulate the distributed functioning of DOC, POC, autotrophic biomass, and the ensuing ecosystem metabolism at the whole network scale. We leveraged on data availability in 12 reaches of a prealpine stream network, the Ybbs River network (Austria), to evaluate the model capability in reproducing local patterns of GPP and ER estimated from the analysis of DO data (Ulseth et al., 2018). We further investigate network-scale emerging properties such as spiraling length and ecosystem efficiency (EE; Newbold et al., 1982;J. Webster & Meyer, 1997).
The mechanistic base of the developed model enables us to investigate how variations in environmental drivers may impact ecosystem functioning and to predict responses of river ecosystem metabolism under changing climatic regimes. Indeed, alterations in either the temperature or the flow regime, or in their mutual interaction, affect OC dynamics and consequently the metabolic regime, shaping the energy pathways, and the EE (Acuña & Tockner, 2010). Warmer temperatures are expected to increase ER to a greater extent compared to GPP (Acuña et al., 2008;Sand-Jensen & Pedersen, 2005), leading to shorter OC spiraling lengths and higher EEs. More intense floods are expected to increase OC transport rates and thus spiraling, to decrease the overall EE and SEGATTO ET AL.

10.1029/2022WR034062
3 of 21 to play a key role in regulating biomass density (Acuña & Tockner, 2010;Raymond et al., 2016;Uehlinger & Naegeli, 1998;Uehlinger et al., 1996), while more pronounced low-flow conditions are predicted to reduce the processing of OC and consequently the metabolic fluxes (Acuña & Tockner, 2010). To test these hypotheses, we investigated how the modeled case study responds to changes in the thermal and hydrological regimes.

Network-Scale Stream Metabolism Model
We propose a spatially explicit, network-scale model of ecosystem metabolism and related OC fluxes (Figure 1). The stream network is discretized into N reaches, that is, the channel segment between two consecutive confluences or between the head of a first-order stream and the next confluence. Each reach constitutes a node of an oriented graph with edges following the flow direction (Bertuzzo et al., 2017;Carraro et al., 2018). The model assumes a stream perspective and thus classifies as autochthonous resources those produced within the stream by autotrophic organisms, which are here assumed to be dominated by benthic communities. All other OC resources, including those produced in the adjacent riparian zone, are considered as allochthonous. Within each node, the model simulates the temporal dynamics of three well-mixed OC pools (DOC [gC], POC [gC], and autotrophic biomass B A [gC]), main biophysical fluxes among them, metabolic fluxes of GPP and ER, terrestrial DOC and POC subsidies, and the advective transport along the flow direction ( Figure 1). The goal is to account for the most relevant processes while maintaining a parsimonious approach and limiting the number of parameters. To this end, we model a single POC pool that includes POC stored in the near surface sediments, POC SED ; POC in The reach-scale ecosystem model describing the temporal evolution of dissolved organic carbon (DOC, light-blue box), particulate organic carbon (POC, yellow box), and autotrophic biomass (B A , green box; Equations 1-3). Organic carbon (OC) mass fluxes originating in the surrounding catchment along with the advective transport from upstream to downstream are represented in orange, while in-stream biophysical processes are shown in black. POC (yellow box) includes all heterotrophic organisms (B H ), and suspended (POC SUS ) and sediment (POC SED ) POC that has entered the system within the simulation horizon. POC BURIED accounts for buried POC and woody debris (gray shaded box) that accumulated mostly outside of the simulation horizon considered and contributes to ecosystem respiration (ER) through the respiratory flux R Basal . Burial (ϕ B ), degradation (L DOC ), loss (L A ), uptake (U DOC and B A ), scouring (S A ), and lysis (L Y ) represent the fluxes from pool to pool, while photosynthesis (gross primary production [GPP]) and all respiratory terms represent reduction of CO 2 to OC and oxidation back to CO 2 , respectively. SEGATTO ET AL.
10.1029/2022WR034062 4 of 21 suspension in the water, column POC SUS ; and heterotrophic biomass (B H ) sustained by this OC pool (Figure 1). POC can be buried into deeper sediments (flux ϕ B ) and contribute to build a long-term storage of OC that contributes to the basal respiration R Basal . Large woody debris is not explicitly modeled as its turn over time could be longer than the time scale of interest of this application (see also discussion below). However, respiration of OC sourced by woody debris can conceptually be included in the basal respiration term, which can be thought of as a term generally accounting for respiration of OC sources that are not explicitly tracked by the model.
For each node i, the model depicted in Figure 1 translates into the following equations: DOC and POC (see Figure 1) can enter the generic i-reach via transport from upstream (ϕ up,DOC,i and ϕ up,POC,i , Equations 1 and 2, note that such fluxes are null for first-order stream reaches); allochthonous lateral inputs in the form of OC leaching from the surrounding hillslopes and riparian zones, and transported by hydrologic lateral flow (ϕ L,DOC,i , Equation  Figure 1 and following Equation 4). The heterotrophic metabolism is also sustained by exudates from autotrophs (e.g., algae) in rivers (Hall & Beaulieu, 2013), flux L A,i in Equation 3. Also in this case, a fraction β produces new heterotrophic biomass and it is thus incorporated into POC ( B A , = A, , Equation 2), while the remaining fraction (1 − β) is spent to perform such a process ( B A , = (1 − ) A, , Figure 1 and following Equation 4). The term R Res,i (Equation 2) encapsulates all other heterotrophic processes resulting in a net oxygen consumption that we did not explicitly include in the model framework, including cellular maintenance metabolism or nitrification, for instance. The outfluxes for the POC pool are represented by downstream advection (ϕ POC,i , Equation 2), burial into the deeper sediments (ϕ B,i , Equation 2), or biophysical lysis and hydrolysis processes (L Y,i , Equation 2), which replenishes the DOC pool. This is, for example, the possible fate of fresh POC from litterfall which can be hydrolyzed and degraded to DOC via physical or biological processes, and then mineralized or incorporated to build new B H mass. Alternatively, litterfall can be transported downstream, as suspended POC, or buried in the streambed. The autotrophic biomass dynamics (Equation 3) is based on the model presented in Segatto et al. (2020): B A,i increases because of the local photosynthetic flux (GPP i ), and it is lost through autotrophic respiration (R A,i ), flow-induced scouring (S A,i ), and by the above-mentioned loss flux L A,i , which encapsulates mortality, detachment, grazing, and exudation (Segatto et al., 2020).
From Equations 1-3, the local ecosystem respiration (ER i ) can be derived as the sum off all respiratory contributions: Finally, the NEP can be calculated as To apply the model, one needs to specify the functional form of each flux. Such a choice should be tailored to the specific case study and depends on the knowledge of the site, drivers expected to control metabolism and crucially on the availability of data for parameter estimation. In the following sections, we first describe the Ybbs River network to then particularize the framework to this specific case study.

River Network Case Study
The Ybbs River network ( Figure 1) drains a 256 km 2 subalpine catchment in Austria (47°48′22.9″ N, 14°57′00.8″ E). The catchment has an average elevation of 937 m.a.s.l. and ranges between 532 and 1,831 m.a.s.l. (Ulseth et al., 2018). The climate is prealpine with an average annual precipitation of approximately 900 mm and an average temperature of approximately 7°C. The area is mainly covered by forests (82%) with minor contributions of alpine meadow (11%) and pastures at lower altitudes (Ceola et al., 2014). The stream network has been delineated from a 10 m resolution digital elevation model (Besemer et al., 2013) and resulted in N = 292 stream reaches. Tree cover density (TCD [-]) has been obtained via the Copernicus Land Monitoring Service. Starting from the 5 of 21 TCD distributed information, we derived for each reach i a proxy of the reach vegetation coverage (TCD r,i [-], i.e., the average TCD over the reach pixels) and watershed vegetation coverage (TCD ws,i [-], i.e., the average TCD over the entire subcatchment contributing to reach i). The former is later used to modulate the spatial distribution of the leaf litter falling directly into the stream reach, while the latter is used as a crude proxy of the organic content of the catchment soils (see e.g., Guo & Gifford, 2002) and is later used to distribute spatially the magnitude of the DOC leaching from soils to stream via lateral discharge.
The catchment has been extensively monitored and studied (Ceola et al., 2014;Segatto et al., 2021;Ulseth et al., 2018), and the relevant data to apply the model are available for the year 2013 (see Table S1 in Supporting Information S1 for a summary of the data set used). Streamflow discharge was measured at the catchment closure at a hourly time step (Q O [m 3 day −1 ]). Stream water temperature (T [°C]) and photosynthetic active radiation (PAR [lux]) were measured in 12 sites and extrapolated at network scale using random forest (Breiman, 2001;Segatto et al., 2021). Daily time series of GPP and ER estimated via the single-station open-channel diel-oxygen method are also available for the same sites (Ulseth et al., 2018) and have been translated to carbon units assuming 1:1 C to O 2 molar ratio (Bott et al., 1978;Burris, 1981;Demars et al., 2016). Data on DOC dynamics (Fasching et al., 2016), benthic biomass (Segatto et al., 2020), and stream geometry (Ceola et al., 2014) were also instrumental to inform the model as detailed in the following.

Model Setup
Each reach i is modeled as a well-mixed reactor with an approximately rectangular cross section of constant width w i (m), length l i (m), time-varying water depth z i (t) (m), and therefore water volume . Given that the Ybbs catchment does not experience strong climatic gradients, discharge in each reach (Q i (t) [m 3 day −1 ]) is assumed proportional to drainage area and has thus been obtained scaling the measured discharge at the outlet (Q O (t)) accordingly to: Segatto et al., 2020), where DA i is the drainage area of reach i and DA O that of the outlet. The geomorphic properties of the Ybbs River network have been extensively studied by Ceola et al. (2014) and we adopted the same exponents to characterize how stream width w i and the Gauckler Strickler's roughness coefficient K s,i (m 1/3 s −1 ) scale as power law functions of the drainage area DA i (Segatto et al., 2021). Water depth z i (t) has been estimated assuming uniform flow conditions via the Manning equation (as in Segatto et al., 2020), which also allows the derivation of the bottom shear stress τ (Pa) as the product between the hydraulic radius r i (m), the specific weight of water γ w (N m −3 ), and the channel slope s i (-) (τ i (t) = γ w r i (t)s i ). Lateral discharge flowing into each reach Q L,i can be derived by continuity as The well-mixed assumption allows relating the advective transport of DOC to its concentration as DOC, = DOC (Equation 5). The DOC flux associated to the discharge coming from upstream reaches thus reads as the first term of the right-hand side (RHS) of Equation 5. Lateral DOC input (ϕ L,DOC,i ) is modeled as the product between lateral discharge Q L,i (t) and its DOC concentration. The latter depends on the specific (i.e., per unit of catchment area) lateral discharge q L,i (third term of the RHS of Equation 5), as widely experimentally observed (see Raymond et al., 2016, for an overview). In particular, Fasching et al. (2016) studied DOC dynamics in the Ybbs catchment and reported strong hydrological controls over patterns of concentration. Accordingly, we adopted the same scaling exponent of ζ = 0.4, while the parameters (a D [gC m −3 ] and b D [gC m −(3+ζ) day ζ ]) have been inferred from measured data in Oberer Seebach, a second-order stream of the Ybbs network (see Table 1 and Fasching et al. [2016]). DOC concentration via lateral flow is further modulated by the average tree cover density in the subcatchment (TCD ws,i ) as discussed above. The removal of DOC (L DOC,i (t)) is parametrized in terms of uptake velocity at 20°, R DOC (m day −1 ; Equation 5), a widely adopted metric in the literature (Mineau et al., 2016). The corresponding removal rate (day −1 ) is defined as the ratio between removal velocity and water depth z i (fourth 6 of 21 term of the RHS of Equation 5). We further account for temperature dependence of the removal process through the exponential Arrhenius model (Jørgensen & Bendoricchio, 2001) with sensitivity θ DOC (-). Lysis of POC into DOC is modeled as a linear process of rate μ Ly (day −1 ; last term of Equation 5).  Note. As not all parameter combinations lead to physically or biologically meaningful model output, we estimated some parameters as a combination of other parameters or variables (fourth column). Range and value (fifth and sixth columns) refer to the parameter combination. If no range is indicated, the parameter value has been fixed based on previous research. Symbols ̄T ,max and ̄L ,max refer to the average (across all reaches) of the maximum (i.e., at maximum local T and PAR) values of the temperature and light limiting functions, respectively. A. stands for autotrophic.

of 21
The advective transport of POC involves the processes of resuspension, transport in the water column, and sedimentation. Assuming an instantaneous equilibrium between suspended and sediment POC, it is possible to express the transport of the bulk POC i (t) pool, which comprises both components, as a function of flow velocity v i and bottom shear stress τ i , which controls the resuspension rate (first two terms of the left-hand side of Equation 6, see Appendix A for the full derivation). High flows increase resuspension rates, which translate into an increased POC flux toward downstream, and vice versa during low-flow conditions. Litterfall (third term of RHS of Equation 6) per unit of streambed area A i (m 2 ) and for a fully canopy covered reach (i.e., for reach tree cover density equal to one, TCD r,i = 1) is modeled as the sum of a constant baseline contribution, which amounts to an annual load of LF b (gC m −2 ; nd y = 365 days to convert annual to daily rates), plus a seasonal contribution which accounts for the increased litterfall during autumn. To mimic the latter, we used a Gaussian function: , where μ LF and σ LF quantify the location and dispersion of the seasonal peak, respectively, and LF y its annual load. A temperature-driven (with sensitivity θ POC [-]) first-order kinetics models the residual respiration R Res,i , with specific rate R POC (day −1 ). A fraction α (-) of the DOC uptake is incorporated into POC as new heterotrophic biomass. POC inputs deriving from autotrophic biomass B A,i are described further below.
Equation 6 tracks the fate of the POC entered and processed during the simulation horizon. In addition to this pool, large woody debris and streambed sediment could store additional OC (generally termed POC Buried in Figure 1) which mostly accumulated before the simulation period. As the simulation for this case study is short (i.e., 1 year long) compared to the time scale of the formation of this pool, we do not specifically model its dynamics (i.e., the initial mass and the flux ϕ B ) but simply the contribution R Basal,i (gC day −1 ) of this pool to ER: where R Basal,O (gC m −2 day −1 ) represents the basal respiration flux per unit of streambed area at the catchment outlet. We assumed the same temperature sensitivity of POC, θ POC . We further adopted a power law function of the drainage area with exponent η (-) to simulate the possible decreasing contribution (η < 0) of this pool from headwaters to larger rivers .

Autotrophic biomass changes in time are expressed as
GPP i is modeled as a function of local active radiation (PAR i (t)) and stream water temperature (T i (t)). μ P (day −1 ) represents the specific photosynthetic rate at 20°C without light limitation (Segatto et al., 2020). Temperature dependence is modulated by the sensitivity θ A (-), while light dependence is modeled through a Michaelis-Menten curve where K PAR (lux) represents the half saturation light (Segatto et al., 2020;Uehlinger et al., 1996). Autotrophic respiration R A,i is controlled by the specific rate μ R,A (day −1 ). For model parsimony, temperature sensitivity of respiration is assumed to be the same of GPP (Segatto et al., 2020). Following Segatto et al. (2020), who successfully applied this model to four reaches within the Ybbs River network, we assume that net growth (GPP i − R A,i ) is modulated by a density dependent term (last factor of the first term of the RHS of Equation 8) which decreases metabolic processes as biomass approaches the carrying capacity D, = D,O(DA ∕DAO) . The latter could change along the network (i.e., it is modeled as a power law function of the drainage area with exponent γ [-]) to possibly mimic the larger carrying capacity of downstream reaches. K D,A,O thus represents the carrying capacity at the outlet. Autotrophic biomass feed the POC pool through loss processes (L A,i (t), including all processes not directly linked to stormflow) and detachment induced by scouring (S A,i (t)). Both have been modeled as linear processes with specific rates μ L,A and μ S (day −1 ), respectively, but scouring is active only when a critical shear stress threshold τ 0,i is exceeded ( The model is solved numerically using a forward Euler scheme with an adaptive time step (longer allowed time step is Δt = 1 hr). To avoid the estimation of initial conditions, we ran a 2-year-long spin-up period by replicating, on an annual basis, the discharge, light, and water temperature time series. Model parameters are summarized in Table 1. Some parameters have been estimated based on empirical data available for the Ybbs catchment or on literature values (Table 1), the remaining 17 parameters have been estimated contrasting simulated and measured time series of metabolic fluxes: GPP (Equation 8) and ER (Equation 4).
Parameter estimation was performed using a two-stage Monte Carlo technique. First, we explored parameters related to autotrophic biomass B A (n = 7, Table 1). Model performance was measured against simulated GPP (as the root mean squared error, RMSE GPP [gC m −2 day −1 ]) in the 12 reaches where metabolic estimates are available. Note that the remaining parameters do not affect the dynamics of the autotrophic biomass B A and hence of GPP. This allowed us to turn off all other processes and reduce considerably the computing time (the most computational demanding processes are the advective transport of both POC and DOC). We thus ran 3 × 10 6 model simulations with parameters randomly extracted from uniform distributions within the range shown in Table 1. Successively, we prescribed autotrophic biomass B A dynamics using the best parameter set found in the first stage, and we explored the remaining n = 10 free parameters against observed ER in the same reaches performing 2 × 10 5 simulations. Following a generalized likelihood uncertainty estimation approach (Beven, 2006), we selected as behavioral the first 5th percentile best performing parameter sets.
Results are presented as spatial and temporal patterns of reach-scale quantities (e.g., areal concentration or fluxes for each stream reach) and network-scale quantities. The latter are computed integrating the reach-scale quantities over all reaches comprising the Ybbs network.

Energy Use, Consumption, and Efficiency
In ecosystem analysis, a useful indicator is the overall efficiency at which ecosystems use available energy (J. Webster & Meyer, 1997). Compared to terrestrial ecosystems, the efficiency of stream ecosystems is more complex to estimate, primarily because downstream transport, and both allochthonous and autochthonous inputs need to be tracked to properly assess OC processing (Fisher & Likens, 1973;Newbold et al., 1982). Our model framework lends itself for this purpose. We calculated, for each reach of the Ybbs River network and on an annual basis (i.e., each flux required for the estimation has been integrated over 1  EE measures the extent to which stream energy inputs are converted to heat and CO 2 within the stream (Fisher & Likens, 1973). Specifically, EE is expressed as the ratio between the amount of OC respired and the sum of both autochthonous and allochthonous OC inputs. In mathematical terms According to our model formulation, we do not consider the contribution of the basal respiration (R Basal ) as it represents the respiration of OC inputs occurred outside of the simulation window and not tracked by the fluxes considered at the denominator. EE ∈ [0, 100%].
Coupling OC mineralization and immobilization dynamics with downstream transport leads to the concept of spiraling (Wallace et al., 1977;J. R. Webster & Patten, 1979), and of spiraling length (Newbold et al., 1982). Spiraling length has two components: turnover length (S T ), which is the average or expected downstream distance traveled by a OC molecule during its residence in the stream, and uptake length (S W ), which is the average distance traveled by an OC molecule in dissolved form in the water column before uptake. Newbold et al. (1982) showed that turnover length S T can be calculated as the ratio between the downstream flux of OC per unit stream width and respiration (cleared from the basal contribution according to our formulation, Equation 10), while DOC uptake length S W can be calculated as the ratio between the flux of DOC per unit of stream width and its uptake rate per unit of streambed area (Equation 11): Spiraling is directly correlated with discharge and it is expected that small streams have shorter turnover lengths (i.e., they are more efficient in using OC) compared to deeper and larger streams where OC is more likely to be transported downstream rather than being mineralized.
Finally, the GPP to ER ratio (GGP/ER) is a measure of heterotrophy (J. Webster & Meyer, 1997) and indicates whether an ecosystem is a net producer or consumer of OC. If GGP/ER > 1, an ecosystem is accumulating OC and it is considered autotrophic; if GGP/ER < 1, the ecosystem is heterotrophic and it is degrading more OC than what is being fixed (Odum & Barrett, 1971). However, it should be noted that in streams GPP/ER depends on import and export and it does not show the extent to which consumers are supported by autochthonous or allochthonous OC (Rosenfeld & Mackay, 1987).

Climate Change Scenarios
We further assessed the impact of changes in temperature and discharge to the simulated metabolic rates of the Ybbs River network. To that end, we ran multiple 1-year-long simulations using the estimated parameters but forcing the system with different thermal and/or hydrological regimes. We explored both independent and joint temperature T and discharge Q variations. We synthetically modified temperature ( * ) and discharge ( * ) hourly time series as displayed in the following equations: We explored scenarios where stream water temperature has been incremented up to T incr = +5°C (Equation 12), whereas we applied an exponential transformation to discharge (with 〈Q i (t)〉 and ⟨ exp ( )⟩ representing the average discharge of reach i under regular and modified scenarios, respectively, Equation 13) where we let vary the exponent Q exp ∈ [1, 2.5] in order to amplify floods and strengthen drought periods. The particular formulation of Equation 13 maintains the same average discharge under the different scenarios to allow for an easier comparison across different runs. Such a choice does not necessarily reflect what is predicted to happen in the future but it allows investigating the isolated effect of changes in allochthonous inputs and OC processing induced solely by more extreme events rather than by an increased or decreased average discharge.

Reach-Scale Organic Carbon and Metabolic Patterns
The model was able to reproduce fairly well the metabolic fluxes in the 12 study sites. Model performance, expressed as RMSE, in reproducing GPP and ER was RMSE GPP = 0.24 (gC m −2 day −1 ) and RMSE ER = 0.28 (gC m −2 day −1 ), respectively (see Figure S4 in Supporting Information S1 for RMSE and normalized RMSE of each site). A summary of the estimated parameters is reported in Table 1. Distribution of the behavioral parameter sets and an assessment of their sensitivity can be found in Text S1 in Supporting Information S1 ( Figures S2-S3 in Supporting Information S1). . The large range of spatial variability of POC is primarily attributable to the low streambed slope of few reaches in the network, which translates into low stream velocity and high accumulation rates.
Results of simulated autotrophic biomass are in agreement with those found in Segatto et al. (2020) in four reaches of the same network. Autotrophic biomass exhibited clear seasonal patterns driven by light and temperature, and generally maintained close-to-equilibrium conditions (i.e., biomass close to carrying capacity, see Figure 2d). The growth season started in late-spring/early-summer where the less limited photosynthetic rate boosted productivity. Biomass started decreasing at the end of the year when climatic conditions became unfavorable. Scouring was activated only during extreme events (see Figure 2d, January 2013 flood event). Maximum autotrophic biomass per unit area was found at the outlet where it was set to be the largest. Average areal concentration was 7.2 gC m −2 (min: 2.1 gC m −2 , max: 10.9 gC m −2 ).

The temporal integration of areal GPP (Equation 8) and ER (Equation 4)
allowed deriving their reach-scale regimes, and by difference the NEP, across the entire Ybbs River network (Figure 3). The model predicted increasing mean daily GPP moving downstream from headwaters (where GPP ranged between 0.04 and ∼0.5 gC m −2 day −1 ) toward the outlet (up to 1.04 gC m −2 day −1 ; Figure 3). On the other hand, mean daily ER was highest in headwaters, especially in the north-and south-west catchments which drain the more forested area and hence have the largest simulated allochthonous inputs. Minimum ER was reached in midsized streams, while ER increased again toward the outlet due to the increased contribution of autochthonous sources: that is, autotrophic respiration and heterotrophic respiration of autotrophic biomass. Overall ER ranged between 0.5 and 1.6 gC m −2 day −1 (Figure 3). The Ybbs River network showed pronounced heterotrophy everywhere (mean daily NEP ranged between −1.6 and 0 gC m −2 day −1 ) except in the main stem, which resulted autotrophic at annual scale (mean daily NEP ranged between 0 and 0.2 gC m −2 day −1 ). In headwater streams, such heterotrophy was mainly driven by basal respiration ( Figure S5 in Supporting Information S1), pointing to a relevant stock of latent OC. Assuming that this respiration flux is mostly sustained by buried POC in the sediment, it is possible to provide a rough estimate of such OC stock. Indeed, if the respiration rate of buried POC ranges between 10 −2 and 10 −3 day −1 (see Catalàn et al., 2022, for an overview of POC respiration rates), the estimated basal respiration implies a stock of buried POC in small streams (DA < 1 km 2 ) ranging between 60 and 1,500 gC m −2 , while the entire network should store on average between 50 and 515 gC m −2 (see Figure S6 in Supporting Information S1).

Network-Scale Organic Carbon and Metabolic Patterns
Spatial integration of reach-scale metabolic fluxes and of OC pools unveiled the ecosystem metabolic regime at the network scale and this across an entire year (Figure 4). GPP peaked in spring ( Figure 4a) and turned the river network autotrophic during such time window (Figure 4c). Throughout the rest of the year, ER dominated (NEP < 0). Network-scale ER resulted partitioned as follows (Figure 4b): 15% DOC respiration (R DOC ), 19% autotrophic respiration (R A ), 24% heterotrophic respiration of autochthonous OC sources ( , 2% of residual POC respiration (R Res ), and 40% of basal respiration (R Basal ), that is, 60% of ER originated from the respiration of OC that entered, or was produced in, the system within the yearly simulation horizon. During the study period, the amount of carbon fixed through photosynthesis equaled 0.23 Gg C (giga-grams of carbon), while respiration, as the total annual ER amounted to 0.31 Gg C. Consequently, NEP at network scale equaled −0.08 Gg C and overall the Ybbs was heterotrophic.
At network scale (Figure 4d), POC was the most abundant OC pool (average, min, and max: 60.1, 1.5, and 221 gC m −2 , respectively), followed by autotrophic biomass (average, min, and max: 15.5, 5, and 23 gC m −2 , respectively). The extreme discharge event of January 2013 mobilized most of the benthic OC. As a result of the fundamental difference in drivers of POC and autotrophic biomass (i.e., temperature, discharge, and resources for the former; temperature, discharge, and light for the latter), recovery pathways were different for the two OC stocks. Autotrophic biomass had, overall, enhanced primary production due to increasing light and temperature with the onset of spring and due to the decreased competition for resources (i.e., biomass far from carrying capacity, see Equation 8). Contrarily, POC recovery, being mostly supported by primary production and allochthonous inputs, was delayed in late summer and autumn, when autotrophic biomass and litterfall were maximum. Expectedly, DOC storage was low compared to the transported fraction and averaged 0.6 gC m −2 (min, max: 0.09, 8 gC m −2 ). top-left), ecosystem respiration (ER; top-right), and net ecosystem production (NEP; bottom-left) as difference between GPP and ER. Plots (bottom-right) show the comparison between the time series of observed (i.e., estimated via the single-station approach from measures of dissolved oxygen) and predicted daily GPP and ER for two representative reach sites (a and b, see location on the maps). The comparison for the remaining sites is reported in Figure S4 in Supporting Information S1. SEGATTO ET AL.
10.1029/2022WR034062 12 of 21 Figure 5 depicts the diverse dynamics of upstream and downstream reaches. Reaches have been clustered and their metabolic fluxes aggregated according to drainage area thresholds. Each bin contains the same number of reaches. GPP increased downstream, while ER decreased downstream to eventually increase toward the outlet as a result of elevated GPP in these reaches. Contribution of autotrophic (R A ) and heterotrophic respiration of the autotrophic biomass ( to ER was proportional to primary production and followed the same pattern, with A typically greater then R A . Respiration of DOC (R DOC ) was maximum in headwaters even though its contribution was noticeable in all reaches. Residual respiration (R Res ) was comparatively low almost everywhere, while the contribution of basal respiration to ER was dominant in headwaters and midsized streams. This scenario depicts a possible longitudinal shift in OC sources contributing to ER, as also predicted by Hotchkiss et al. (2015): from headwaters where respiration is sustained by the surrounding hyporheic and riparian zone to larger streams where respiration is benthic dominated. Figure 6 shows the EE metrics of the Ybbs's River reaches as a function of their position along the network, here indicated by catchment size, and for which we established and explored the emerging scaling relationships. Spiraling length (Figure 6a) increased downstream, which was expected given that headwaters are typically more reactive (Battin et al., 2008). DOC uptake length (S W , average, min, and max: 103, 5, and 730 km, respectively) was lower than spiraling length (S T , average, min, and max: 64, 8.9, and 234 km, respectively) in small streams with drainage area less than 1.6 km 2 ; this indicates that in downstream reaches, DOC was more rapidly transported than stored into biomass. The ratio GPP/ER (Figure 6b) increased sharply for catchment size up to 4.9 km 2 , beyond which it kept increasing but at a slower rate, eventually becoming positive toward the outlet. GPP/ER equaled one in the main stem at DA = 103 km 2 and was maximum at catchment outlet (GPP/ ER average, min, and max: 0.34, 0.02, 1.24, respectively). EE (Figure 6c), on the other hand, did not change significantly with catchment size below 1.6 km 2 (interestingly equal to the intersection point between uptake and turnover lengths), but beyond this threshold, EE significantly decreased as a power law function of catchment size. The simulated EE values depict reaches in which OC that enters during 1 year is poorly processed within the same time window (EE average, min, and max: 1.5%, 0.1%, 10%, respectively). OC inputs are thus preferentially exported from the reach and/or stored at annual scale. We note that if one includes basal respiration in the EE calculation, that is comparing the annual inputs with the total amount of OC respired during the same year (as customarily done in literature, see e.g., Fisher & Likens, 1973), EE increases on average everywhere in the network (form 1.5% to ∼6.5% at network scale), but especially in low-order reaches where it can locally peak at 50% (see Figure S10 in Supporting Information S1). Figure 7 shows the variability of the network-scale metabolism induced by simulating warming stream water temperature and an increased streamflow flashiness (see Figure S7 in Supporting Information S1 for the induced effect on the network-scale OC stocks and Figure S8 in Supporting Information S1 for the averaged network effect on metabolic fluxes). Temperature affects biological rates at different magnitudes according to their sensitivities. As a result of the faster degradation and primary production rates, temperature increments led to higher metabolic fluxes throughout the entire network (Figures 7a and 7b). GPP at maximum alteration (T incr = 5°C) increased by 60% (0.36 Gg C) compared to the baseline scenario, while ER increased by 80% (0.56 Gg C). Therefore, the  (Figure 7c) was an increased OC consumption and potentially increased CO 2 outgassing (cumulative NEP at maximum temperature alteration equaled −0.2 Gg C, i.e., 150% increase compared to baseline scenario).

Climate Change Scenarios
The effect of an altered hydrological regime is more complex as it nonlinearly affects advective transport, allochthonous inputs, mobilization, and scouring rates. Autotrophic biomass and POC were progressively slower in recovering from the progressively stronger flood event of January ( Figure S7 in Supporting Information S1) and this negatively affected GPP in that period (Figure 7d). On the other hand, a flashier hydrologic regime led to two new localized scouring events in June and July, when, at maximum disturbance (Q exp = 2.5), GPP peaked at 2.14 and 1.83 g C m −2 day −1 , respectively. This effect is a direct consequence of the autotrophic biomass being pushed below its carrying capacity (Equation 8) thus leading to an increased rate of growth, a well-known feature of recovering ecosystems (Dzubakova et al., 2018;Odum, 1969). The annual cumulative GPP equaled 0.23 Gg C (0% relative variation, see Figure  S8 in Supporting Information S1). DOC and advected POC followed discharge dynamics (see Figure S7 in Supporting Information S1), increasing during high flow events and diminishing during droughts. This pattern translated into lower resource availability (i.e., higher transport rates during high flows and less DOC input during droughts) which led to a slight decrease in ER (Figure 7e; cumulative ER at maximum discharge alteration: 0.25 Gg C, −20% relative variation), which however was limited in magnitude as ER was sustained by the basal respiration which is unaffected by discharge. Consequently, the overall effect on NEP (Figure 7f) was toward a more autotrophic network caused by the new positive NEP windows induced by scouring events is spring and summer when conditions for primary production were ideal (cumulative NEP at maximum alteration: −0.02 Gg C and 75% relative decrease).
When combining both temperature and streamflow variations in the same simulation, the result was intermediate. GPP evolved in the same way as it did when discharge was altered ( Figure 7g); however, with the difference that maximum rates were higher due to the higher temperature (at maximum variation, annual cumulative GPP = 0.36 Gg C, 56% relative increase). The discharge effect on ER was overridden by that of increased temperature (Figure 7h), mostly because higher temperatures affect also basal respiration (at maximum alteration, annual cumulative ER = 0.46 Gg C, 48% relative increase). The strong tendency of the Ybbs River network to become more heterotrophic because of higher temperatures was attenuated by the increased GPP (Figure 7i) especially immediately after new flow-induced disturbances (at maximum variation, annual cumulative NEP = −0.1 Gg C, 25% relative decrease).
Network-scale energy processing changed predictably as a consequence of our simulated climate change scenarios ( Figure 8): spiraling lengths, both in terms of uptake and turnover lengths, decreased almost linearly with increasing temperature (due to higher DOC uptake and respiration rates) and increased with increased flash iness (due to higher advective fluxes), resulting to a net increase when simultaneous variations were tested (Figures 8a  and 8b). Temperature effect on GPP/ER, as stated above, was in favor of ER (i.e., the ratio GPP/ER decreased), while discharge had a lower impact on ER compared to GPP which, after the emergence of new scouring events, saw increased production rates. As a result, GPP/ER first decreased for low discharge variations but then increased again for high variations (Figure 8c). The combined effect was toward smaller GPP/ER but mitigated by the hydrological regime. Finally, the Ybbs River network became more efficient in using available energy when temperature increased, due to the stronger impact on ER than on GPP, yet less efficient when advective transport was boosted, and intermediate but still less efficient when temperature and discharge were simultaneously changed (Figure 8d).

Discussion
Both biomass and metabolic processes that characterize any ecosystem are constrained by the total amount of energy that is either fixed within or delivered across its boundaries . Since the introduction of the River Continuum Concept (Vannote et al., 1980), our understanding of the spatial and temporal patterns of the metabolic processes and fluxes in streams has steadily advanced. However, we have not been able to assess metabolic regimes at a scale large enough to integrate the effects of hydrological connectivity and environmental drivers on stream ecosystems. Clearly, this is at the scale of entire stream networks. In this context, there is also rising interest on how climate change (i.e., elevated temperature and more extreme hydrological regimes) will affect freshwater ecosystems, energy pathways, and ultimately regional and global carbon budgets.

The Ybbs Network Model
The proposed model provides a simple meta-ecosystem framework for the estimation and process-based simulation of metabolic regimes over an entire stream network. We have shown how interactions among different OC pools (e.g., DOC, POC), autotrophic biomass, and environmental variables (e.g., temperature, PAR, discharge, and geomorphology) were instrumental to simulate consistent patterns of GPP, ER, and NEP across scales in a real-world stream network. As for any modeling exercise, we stress that results depend, albeit with different sensitivity, on the process parameters either estimated or assumed, on the model structure (i.e., the functional form of each flux), and on the observational data sets used for parameter estimation. While this should caution against overconfidence in exact values and predictions derived from our model results, it also provides new motivation, direction, and context for future improvement and development, as discussed in this section.
The value of the parameters controlling the functional forms of the DOC concentration entering via lateral flow (a D , b D , and ζ) was assumed based on a scaling relationship derived for a reach of the same network of our case study (Fasching et al., 2016), while those controlling the annual litterfall load and its timing (LF b , LF y , μ LF , and σ LF ) according to published data in Bretschko (1990) and Neumann et al. (2018). The spatial variability of these two fluxes has been modulated by remotely sensed information on TCD at reach and subcatchment scales. We argue that the current framework could potentially infer more accurate patterns if coupled with an explicit hydrological model and channel flow-routing scheme (e.g., Payn et al., 2017), enabling the direct estimation of lateral discharge, water stage, flow velocity, and other relevant hydrological fluxes (e.g., snowmelt) at any reach of the network. The coupling with a soil submodule describing major terrestrial processes (e.g., soil carbon leaching) would also enhance the accuracy (Grandi & Bertuzzo, 2022). This improvement would allow relaxing most of the assumptions detailed in Section 2 and more properly simulating the spatial variability in allochthonous inputs. Moreover, the application presented here uses the stream reach as the basal spatial unit. However, the same framework can directly be applied at a finer spatial resolution (e.g., a regular one-dimensional grid) if useful information to characterize within-reach heterogeneity is available.
We do acknowledge that this specific application failed to reconstruct some localized events of ecosystem metabolism. For instance, model results consistently underestimated the boosted GPP observed in some reaches during spring (see e.g., sites 2, 6, and 12 in Figure S4 in Supporting Information S1) or the spiked ER during the same period or during high flow events in fall (see e.g., sites 2, 3, and 6 in Figure S4 in Supporting Information S1). This mismatch could be related to incorrect characterization of environmental variables or vegetation cover, to site-specific differences and heterogeneity in community composition, or to the flushing of nutrients during snowmelt (see e.g., the analysis of Ulseth et al., 2018, on the same catchment). Also the analysis performed by Segatto et al. (2021) using a Random Forest algorithm on the Figure 6. Scaling of ecosystem efficiency metrics as a function of drainage area. Plots are in log-log scale and each dot represents one stream reach. Scaling relationships have been derived using piece-wise regression lines where breakpoints have been selected in order to minimize the overall root mean squared error (RMSE) among all possible combinations (including no breakpoints). (a) The spatial variability of both uptake length (S W , pink) and turnover length (S T , green) averaged over the entire simulation period along with their regression lines (no breakpoints individuated, S W = 33.02 × DA 0.58 , p < 0.001, R 2 = 0.99; S T = 35.17 × DA 0.39 , p < 0.001, R 2 = 0.97). Yellow dot highlights the intersection point between the regression lines. (b) Annual averaged gross primary production (GPP) to ecosystem respiration (ER) ratios (blue) and the corresponding piece-wise regression lines (breakpoint located at DA = 4.9 km 2 ; lbp = left side; rbp = right site of the break point: 10.1029/2022WR034062 15 of 21 same data set revealed a seasonal contribution to GPP and ER that could not be explained by the variation of the other environmental variables. These dynamics are currently not accounted for in our framework, but could potentially be accommodated, at least in part (e.g., for snowmelt), by developing the coupled model described above. Figure 7. Network-scale metabolic response to alterations in thermal and hydrological regimes. Temperature has been increased up to +5°C, whereas the exponent of discharge transformation has been varied between 1 and 2.5. Each colored line represents a trajectory that the network-scale metabolic flux performs under a specific alteration (see colorbars). Warming effects to network-scale gross primary production (GPP), ecosystem respiration (ER), and net ecosystem production (NEP) are displayed in (a)-(c); variability induced by discharge alteration is shown in (d)-(f). (g-i) The combined effect of both temperature and discharge alterations. Refer to Figure S7 in Supporting Information S1 for the induced variability to dissolved organic carbon (DOC), particulate organic carbon (POC), and B A . The application presented here was made possible thanks to a synergy of different monitoring campaigns carried out in the Ybbs River network (see Table S1 in Supporting Information S1 for a summary of the data set used), and this poses challenges to a straightforward transferability of the model to other case studies. The required suite of data set depends on the dominant OC sources and pathways in the studied network and it could be less demanding than the one presented here. For instance, in this application a lot of effort was dedicated to characterize DOC input through the analysis of the available high-frequency (subdaily) DOC measurements, information typically not available in standard water quality monitoring programs. However, the estimated contribution of DOC to the network metabolism was minor and with low temporal and spatial heterogeneity. This suggests that coarser information on this OC pool could have sufficed. So it is advisable to first identify the dominant OC sources in order to guide the experimental sampling effort. To this end, a preliminary analysis could be carried out following the method proposed by Bertuzzo et al. (2022), which allows partitioning the contribution of different OC sources to ER using DO time series. Regarding other environmental variables, the acquisition of data in any catchment is progressively facilitated by recent methodological advances which combine satellite information with modeling, see for instance Revel et al. (2021) and Hossain et al. (2022) for streamflow discharge, Savoy et al. (2021) for light availability at streambed, and Bertuzzo et al. (2022) for litterfall input derived from Leaf Area Index remote observation. Finally, crucial to the applicability of the presented model are long (annual or multiannual), high-frequency time series of DO, which are increasingly available thanks to the recent advent of cheap sensor technology . We stress the importance of designing a nested sampling scheme to be able to capture the scaling of relevant variables (see e.g., exponents γ and η, Table 1) with catchments size.

Multiscale Stream Ecosystem Metabolism and Efficiency
The model predicted increasing mean daily GPP moving downstream from the headwaters toward the outlet, where wider channels, with reduced canopy cover and highest carrying capacity promoted GPP and autotrophic biomass (Figure 3). On the other hand, highest mean daily ER was simulated in the headwaters, where the allochthonous inputs peaked, and basal respiration (R Basal ) was maximum. Minimum ER was observed in midsized streams, while the outlet resulted more heterotrophic due to the increased contribution of the autotrophic respiration (R A ) and of the respiration of autotrophic byproducts ( . This general pattern is in line with earlier predictions on stream ecosystem metabolism (Fisher & Likens, 1973;Hotchkiss et al., 2015;Vannote et al., 1980) and studies compiling data from a wide suite of stream ecosystems (Battin et al., 2008;Bernhardt et al., 2018;Bertuzzo et al., 2022;Diamond et al., 2021;Hoellein et al., 2013).
Our findings confirm previous results (Segatto et al., 2021) and show that the metabolism of the Ybbs River network is heterotrophic at the annual scale. The network-scale GPP regime (Figure 4) was mostly driven by the downstream reaches and showed a conspicuous peak during spring, low activity during winter, and intermediate production during the rest of the year. Our result extrapolated for a real case study resembles that obtained by Koenig et al. (2019) using theoretical optimal channel networks under the stochastic assignment, that is, the scenario proposed to be more representative of real rivers by the authors. The simulated contribution to the network-scale ER was more evenly distributed along the network (Figure 4). However, this more regular pattern is achieved thanks to a clear shift in the OC sources sustaining respiration ( Figure 5). Our annual estimates of GPP and ER at network scale are also closely bracketed by those reported from other studies covering various biomes (e.g., Battin et al., 2008;Hoellein et al., 2013;Rodríguez-Castillo et al., 2019). Moreover, converting these fluxes back to O 2 units (using a 1:1 M ratio), our estimates of GPP, ER, and NEP equal 0.61, −0.83, and −0.22 Gg O 2 , respectively. These figures are in very good agreement with the results shown in a recent machine learning application to the same river network (Segatto et al., 2021), where the authors reported for the same year 0.63, −0.83, and −0.20 Gg O 2 , respectively. We are therefore confident that our framework represents a valuable tool to simulate annual metabolic regimes in every reach of an entire stream network. Moreover, the mechanistic nature of the model allows predictions on the ecosystem responses to climate or land use changes, or to other anthropogenic disturbances.
An important novel aspect of the model is that its structure allowed us also to infer the contribution of different OC sources or pathways to the total ER (Figure 4) approach involving machine learning and scaling of allochthonous resources along the network. On the other hand, we could not trace back the source of the basal respiration that accounted for 40% of ER. This contribution is tentatively attributed to the respiration of OC buried in the streambed sediments, which had mostly built up before the simulated period. Due to the limited duration of the simulation (i.e., 1 year), we could not model the dynamics of this OC stock which changes over longer time scales, but simply accounted for its contribution to ER. This limitation could potentially be overcome with longer data sets and simulations which could unlock also the estimation of burying rates and the fluctuations of this stock. To assess the plausibility of this contribution, we assumed a wide range of possible respiration rates and estimated the stock of buried OC required to sustain the estimated respiration ( Figure S6 in Supporting Information S1). Results are in agreement with other studies assessing benthic and buried POC (Magana & Bretschko, 2003;Smock, 1990;J. Webster et al., 1995).

Climate Change Effects on Ecosystem Functioning
We investigated how variations in the environmental forcings could impact ecosystem functioning and explored in which ways our simulated case study may react if forced with different thermal and hydrological regimes. Sand-Jensen and Pedersen (2005) showed that 4°C-5°C rise in air temperature will increase oxygen consumption rates at ambient temperature by 30%-130%, while Acuña et al. (2008) shown that a 2.5°C air temperature increase will increase ER by an average of 20%. Our model also predicted that at maximum stream water temperature increase (i.e., T incr = 5°C), network-scale ER will increase by 80% (53%-180%), while at a T incr of +2.5°C, an increase of ∼31% (24%-57%) is expected. In line with the estimated contributions to the total ER, a large share of this enhanced respiration is due to basal respiration. However, it should be noted that higher respiration could lead, in the long term, at a depletion of the OC pool buried in the sediment. Also in this case, modeling explicitly the dynamics of this OC stock could better constrain these projections. Similarly to predictions of Yvon-Durocher et al. (2010), who combined experimental field data with theoretical predictions derived from the metabolic theory of ecology, our simulated temperature effect on GPP resulted less marked, depicting the Ybbs River network with higher levels of heterotrophy and potentially higher CO 2 emission with warming temperature (Figure 7 and Figure S9 in Supporting Information S1). Higher metabolic rates with warmer temperature translated, at network scale, in shorter spiraling lengths and increased overall EE, confirming our initial hypothesis (Figure 8). Such predictions critically depend on the estimated values of the temperature sensitivity parameters θs. Best model performance was achieved with θ A = 1.103, θ POC = 1.135. When expressed in terms of activation energy, these values correspond roughly to 0.66 and 0.87 eV K −1 , which are in line with, or slightly higher than those reported earlier (see e.g., Acuña & Tockner, 2010;Acuña et al., 2008;Gillooly et al., 2001). It should also be stressed that in our warming scenarios, lowland reaches, which host the large share of the total streambed area, experienced temperatures higher than those observed anywhere in the data set used for parameter estimation. Extrapolation outside the observed range should always be interpreted with caution.
The response to altered hydrological regimes was more subtle and less intuitive. As highlighted by Miao et al. (2009) and Acuña and Tockner (2010), models should take into account the entire sequence of hydrologic events and their prolonged effect on the living biota in order to assess complex responses to alterations of the regime type. We have shown how more extreme hydrological regimes had more pronounced effects on GPP and particularly on its recovery rates after scouring (Figure 7). We observed that the effects of spates depend on the biomass state, that is, biomass far from or close to its carrying capacity. In the former case, disturbances abated production rates, while in the latter they created new room for autotrophic development with a consequent boost in growth rates and GPP. The effect of a modified flow regime on ER was buffered by basal respiration which was unaffected by the flow, but it was still observable and mostly driven by changes in autotrophic respiration and by the reduced processing of OC induced by the pronounced low-flow conditions limiting OC transport. At network scale, flow alterations directed the Ybbs River network toward more autotrophic states due to the appearance of new positive NEP windows induced by scouring (Figure 7 and Figure S9 in Supporting Information S1), higher transport rates, longer spiraling lengths, and lower EE (Figure 8). The antisynergistic effect of temperature and discharge posed their mixed-effect in between the single end-members. The emergence of more scouring events allowed for higher GPP rates due to simultaneous higher temperatures, while ER was primarily driven by temperature, increasing with warming (Figure 7 and Figure S8 in Supporting Information S1).
It should be stressed that we have simulated only direct impacts of climate variations on the Ybbs River network, without considering biological changes. In fact, shifts of thermal and hydrologic regimes could differently affect ecological communities and their emergent response. Regarding climatic variations, other limitation are identifiable. We did not include, for example, variability in temperature at the annual (Sand-Jensen & Pedersen, 2005) or daily (Dang et al., 2009) scales. The tested hydrological regimes allowed to increase the frequency and magnitude of floods and droughts as typically predicted (Acuña & Tockner, 2010;Pletterbauer et al., 2018), but Equation 13, along with our hypothesis of nearly rectangular cross section, does not permit to assess variations in network connectivity induced by the complete disappearance of flow in some reaches or to account for floodplain interactions (Catalàn et al., 2022). While these caveats and limitations highlight the uncertainty related to this kind of projections, we argue that a relevant result of this study is unveiling the complex set of interactions between biota and the environmental matrix, and its control over the metabolic regime and its possible response to climate changes. A perturbation at one point in time can have lasting effects on both biogeochemical fluxes and OC pools. Alterations of the metabolic balance in headwaters change the flux of material and energy flowing downstream and could be counterbalanced or enhanced by the response of lowland reaches. Our analysis highlights how these processes can be appreciated only going beyond the traditional field of view of stream ecosystem research-the reach, and looking at the system from a network perspective. We stress that the present study was made possible by a data set comprising time series of metabolic metrics at several stations nested within the same network. The recent development of sensor technology has revolutionized our capabilities of measuring metabolism in rivers , and we predict that other similar data sets will soon be available. We argue that the modeling tools used to analyze and interpret this data deluge should keep the pace and that the framework presented herein is an effort in that direction.

Appendix A: POC Transport Dynamics
This section derives the simplified POC transport equation used in the main text. Let consider a well-mixed stream reach of volume V = w × l × z (m 3 ), where w, l, and z (m) are, respectively, stream width, length, and stream water depth, and in which POC is suspended in water (POC SUS [gC]) with concentration [POC] SUS = POCSUS∕ = POCSUS∕( × × ) (gC m −3 ). Let denote with POC SED (gC), the mass of sediment POC and [POC] SED = POCSED∕( × ) (gC m −2 ) its aerial concentration. Let POC (gC) be the total amount of POC stored in the system: POC = POC SUS + POC SED . Suspended POC is subjected to advective transport (ϕ POC ) and sedimentation with specific rate μ sed (day −1 ; Figure 1), while POC SED can be mobilized, that is, resuspended at a rate μ res (day −1 ), before being advected. ϕ POC can be expressed as where v = Q/(w × z) (m day −1 ) represents the average stream water velocity.
Let now assume that the dynamics of the subsystem POC SUS -POC SED , ruled by sedimentation and resuspension processes, are fast enough to guarantee an almost instantaneous equilibrium (separation of fast and slow dynamics): where v sed (m day −1 ) is the sedimentation velocity. Typically, v sed ≫ μ res (Acuña & Tockner, 2010) so that at any time most of the POC is stored in the sediment rather than suspended, we can thus approximate POC SED ≈ POC. where v mov = μ res /μ sed (-) is a dimensionless parameter describing the tendency of POC to be mobilized and transported downstream. This simple approximation allows modeling solely the total mass POC instead of both POC SED and POC SUS and their interactions (sedimentation, resuspension, etc. in Equation 6 is derived modulating v mov according to a power law function of the bottom shear stress in order to account for higher mobilization rates during high flow events.

Appendix B: Autotrophic Respiration Fraction
To estimate the fraction β of the loss term of the autotrophic biomass L A which contributes to the POC build up (i.e., it is not directly respired), we combined the model output for the autotrophic biomass compartment (which calibration is independent of β, see Section 2.3), with the method introduced by Hall and Beaulieu (2013) to estimate the fraction of GPP which respired by autotrophs (flux R A in the proposed model) and their closely associated heterotrophic bacteria (flux B A ). This fraction is termed Autotrophic respiration fraction: ∕GPP and it can be estimated from the slope of the 0.9 quantile regression between daily estimates of ER and GPP (Hall & Beaulieu, 2013). We performed this analysis for the 12 experimental reaches ( Figure S1 in Supporting Information S1), which led to an average AR f of 0.59, in line with the range estimated in Hall and Beaulieu (2013). From the calibrated model (autotrophic compartment), we computed the average ratio of R A /GPP which yielded 0.26. It does follow that the heterotrophic contribution on average must be equal to 0.59-0.26 = 0.33. To achieve this value, we derived, from the modeled autotrophic compartment, that β must be equal to 0.55.

Data Availability Statement
No new data were generated for this study. Data used come from a previous publication (Segatto et al., 2021) and are available at the following public repository https://doi.org/10.5281/zenodo.3972937.