Cross‐Scale Interactions Dictate Regional Lake Carbon Flux and Productivity Response to Future Climate

Lakes support globally important food webs through algal productivity and contribute significantly to the global carbon cycle. However, predictions of how broad‐scale lake carbon flux and productivity may respond to future climate are extremely limited. Here, we used an integrated modeling framework to project changes in lake‐specific and regional primary productivity and carbon fluxes under 21st century climate for thousands of lakes. We observed high uncertainty in whether lakes collectively were to increase or decrease lake CO2 emissions and carbon burial in our modeled region owing to divergence in projected regional water balance among climate models. Variation in projected air temperature influenced projected changes in lake primary productivity (but not CO2 emissions or carbon burial) as warmer air temperatures decreased productivity through reduced lake water volume. Cross‐scale interactions between regional drivers and local characteristics dictated the magnitude and direction of lake‐specific carbon flux and productivity responses to future climate.


Introduction
Inland lakes are key elements of both global freshwater ecosystems and the global carbon cycle. For example, lake primary productivity supports globally important fisheries that are the main protein source for many human populations (Lynch & MacMillan, 2017) and constitute a $25 billion per year recreational fishing industry in the United States alone (U.S. DOI et al., 2011). Collectively, lake carbon (C) fluxes also constitute a significant portion of the global C cycle, as current estimates of C exported annually from terrestrial ecosystems to inland waters are on par with global land net ecosystem production (NEP; Cole et al., 2013;Drake et al., 2017;Randerson et al., 2002;Tranvik et al., 2009).
Climate change has already impacted lakes and reservoirs in heterogeneous ways (McCullough et al., 2019) as cross-scale interactions (Soranno et al., 2014) between climate and lake-specific characteristics have caused some lake water temperatures to increase faster than air temperature and other lake water temperatures to cool despite increased air temperature (O'Reilly et al., 2015;Rose et al., 2016). These lake-specific changes in water temperature impact organismal metabolic rates, with heterotrophy more sensitive to changes in temperature than autotrophy Yvon-Durocher et al., 2010). Watershed inputs of C and nutrients are also likely to change with predicted changes in precipitation. Extreme precipitation events, which have been increasing in intensity in areas like the Midwest U.S. (Min et al., 2011), account for a significant fraction of C and nutrient loads to lakes Zwart et al., 2017), which impact lake productivity (Zwart et al., 2017), CO 2 emissions from the lake surface (Vachon & del Giorgio, 2014), and the percent of terrestrial C removed within lakes due to changes in lake water residence time (WRT) and metabolism (Evans et al., 2017;Jones et al., 2018).
Several single-lake models have been run with future climate scenarios, documenting climate change effects on individual lakes (Kiuru et al., 2018;Tan et al., 2018). However, given the cross-scale interactions between regional climate drivers and local characteristics, integrated modeling of many lakes at regional scales is needed to inform whether lakes collectively might increase or decrease CO 2 emissions and carbon burial and how primary productivity may change. Here, as a regional-scale case study, we use an integrated modeling framework of lake hydrology and biogeochemistry Zwart et al., 2018) to project changes in lake-specific and regional primary productivity and C fluxes under 21st century climate change simulations for thousands of lakes in Northern Wisconsin and Michigan, USA. Our integrated model was driven by statistically downscaled daily temperature and precipitation from a set of six representative global climate models (GCMs) from the CMIP5 archive, all using the highest emissions scenario (RCP8.5), during two separate time periods within the 21st century (2050s and 2080s), resulting in a total of twelve regional simulations representing the future climate (Byun & Hamlet, 2018).

Modeling Overview and Domain
Our study leveraged newly published tools developed to simulate detailed hydrological and biogeochemical fluxes for thousands of lakes and reservoirs over large spatiotemporal scales. The lake hydrology model  utilized a computationally efficient integrated surface water (SW) and groundwater (GW) modeling framework that informed a lake water budget model incorporating daily hydrologic inputs and exports from individual lakes within the modeling domain. The lake biogeochemical model was informed by the hydrologic information and was built upon a simple lake heat budget, constituent loading, and lake biogeochemical model to track carbon storage and processing for all lakes at a daily timescale within the modeling domain ; supporting information Figure S1).
For this current study, we modeled the same lake-rich region analyzed previously in Hanson et al. (2018) and Zwart et al. (2018), the Northern Highlands Lake District (NHLD) located in a forested and pristine portion of northern Wisconsin and Michigan, USA ( Figure 1j). Specifically, we modeled 3,675 lakes and reservoirs within an area of about 6,000 km 2 . Given that these models reasonably estimated hydrological and biogeochemical variables historically  for this region Zwart et al., 2018), we projected the models into the future to investigate the hydrological and biogeochemical response of NHLD lakes under projected future climate scenario periods that are representative of the years 2041-2100.

Meteorological Forcing Data and Selected Climate Scenarios
We forced our integrated model with statistically downscaled daily temperature and precipitation data using the Hybrid Delta (Hamlet et al., 2013;Tohver et al., 2014;HD) approach applied within the Midwest and Great Lakes region of the United States and Canada by Byun and Hamlet (2018). We used six GCM outputs (supporting information Table S1) that (a) performed well across the Midwest and Great Lakes region and (b) represented the range of changes in temperature and precipitation from a larger ensemble of 31 GCMs. Byun and Hamlet (2018) provide additional details on these methods and analyzed 10 GCM ensembles for three time periods and two emissions scenarios (a total of 60 scenarios). Here, we focus our attention on a subset of six representative scenarios spanning the full range of results for the 2050s and 2080s time periods for the RCP8.5 (highest) greenhouse gas concentration scenario. We use the highest greenhouse gas concentration scenario because this greenhouse gas concentration scenario compares well with current emissions trajectories (Hayhoe et al., 2017).

Coupled Surface Water and Groundwater Hydrology Model
For a complete, detailed description of the hydrologic modeling approach, see Hanson et al. (2018). The hydrology model utilized an integrated SW/GW modeling framework that combined the Variable Infiltration Capacity (VIC) land surface hydrology model (Liang et al., 1994), the GFLOW analytic element method groundwater model (Haitjema, 1995), and a lake water budget model. Our modeling framework simulated the daily changes of water storage within each individual lake due to water flux inputs and exports that are informed by the SW and GW models. As VIC is a translator of climate information to water information, the only information that changed within our hydrologic modeling framework was the meteorological forcing data input to VIC with appropriate changes in temperature and precipitation created by Byun and Hamlet (2018). Average monthly air temperature (a) and monthly sum of precipitation (c), evapotranspiration (e), and streamflow (g). Mean regional values across all six GCMs for future climate drivers are shown in the solid lines for the 2050s (orange) and 2080s (red), with the minimum and maximum bounds shown by respective shaded coloring. The historic mean regional monthly driver values are shown in dashed blue line. The average annual difference between the future scenarios and historic model runs for temperature (b) and cumulative sum between the future scenarios and historic model runs for precipitation (d), evapotranspiration (f), and streamflow (h) are shown to the right of each panel for the 2050s (orange) and 2080s (red). Projected future changes in ice cover and snow depth for a single characteristic lake  are show in panel i, where the densities of ice on (points) and ice off (crosses) for the three separate time periods are shown in their respective colored shading along with the longest, average, and shortest ice duration for each time period. Average snow depth by day of year (gray shading) is also shown in panel i with minimum, mean, and maximum snow depth across the six GCMs are overlaid on top of each other, while the historic model run only has one realization of average snow depth by day of year. The range of seasonal maximum snowpack using the six different GCM projections is shown for the future scenarios, while only one maximum snow depth is shown for the historic model run. Panel j shows the study region with the modeled lakes shown as blue polygons.
Model parameterization and simulation methods for the hydrology modeling are identical to those used in Hanson et al. (2018) for the future climate scenarios presented here. The future climate simulations were initialized in the same manner as the historical (i.e., lakes start at the same lake stage elevation) and were performed over water year (WY) 1980(1 October 1979to 30 September 2013. Historical simulations in Hanson et al. (2018) ran through WY 2015, but the Byun and Hamlet (2018) HD data sets used for this study only extended through WY 2013. All lake basin characteristics and morphometry parameterization were the same as detailed in Hanson et al. (2018), using modifications for a subset of lakes (n = 7) described in the discussion section of that work.

Lake Heat Budget, Constituent Load, and Biogeochemical Models
Here, we briefly describe the lake models used for this study, but for more details, please see Zwart et al. (2018). The lake heat budget, constituent load, and biogeochemical models were driven by the hydrologic output from the coupled surface water and groundwater model and parameterized exactly as described in Zwart et al. (2018) with altered driving data from the climate change scenarios as described above. Our lake heat budget model used mass balance equations (Lenters et al., 2005) and functions from the R package LakeMetabolizer  to model lake water temperature during the open-water period as in Zwart et al. (2018). The constituent loading model uses land cover, long-term data, and simulated water budgets from our coupled hydrologic model to estimate inflowing dissolved organic carbon (DOC), dissolved inorganic carbon (DIC), terrestrial particulate organic matter, and dissolved inorganic phosphorus to each lake Zwart et al. (2018). For this analysis, we focused on the impact of climate change scenarios on lake C fluxes and productivity. Since we did not evaluate the impact of land cover change on lake C fluxes and productivity, we kept inflowing constituent concentrations for the future scenarios the same as they were for historic scenarios. The lake biogeochemical model takes watershed inputs of constituent loads along with meteorological driver data and translates the inputs into daily estimates of CO 2 emissions from the lake surface, C burial in the lake sediments, C export through groundwater and surface water pathways, and primary productivity . CO 2 exchange with the atmosphere was a function of the concentration gradient between the surface water and the atmosphere and modeled piston velocity. We projected future daily atmospheric CO 2 concentrations by using projected average annual CO 2 concentrations from the RCP8.5 scenario (Riahi et al., 2011), adjusted by intra-annual variation from the Keeling curve time series (Keeling et al., 2001). This produced intra-annual variation in CO 2 concentration that was representative of the Northern Hemisphere and future CO 2 concentrations that reached~950 ppm by the end of the 21st century.

Model Runs, Aggregation, and Data Availability
The lake biogeochemical models were run in parallel using a Docker container on the HTCondor distributed computing cluster (Erickson et al., 2018), with a total compute time for all 12 scenarios for 3,675 lakes lasting about 5 hr using 1,400 cores. Model results were aggregated from daily output to annual means for each lake and output variable, including CO 2 emissions, carbon burial, primary productivity, and morphological and physiochemical characteristics. Regional C flux and productivity response variables were computed as the sum of lake-specific average annual fluxes. The aggregated model output used for all analyses and figure generation, the process-based model code, and analyses code are available on v1.0 of our GitHub repository (https://github.com/jzwart/NHLD_climate_change; v1.0 release: https://doi.org/10.5281/zenodo.3352291). Full, daily time step model output (>50 Gb) is available from the data release described by Zwart et al. (2019).

Statistical Analyses
We regressed regional C flux and productivity response variables against the cumulative annual difference in precipitation and evapotranspiration. Long-term average precipitation minus evapotranspiration and streamflow are highly significantly correlated across all climate scenarios (p < 0.001, R 2 > 0.99), and streamflow has been shown to impact lake C fluxes and whole-lake metabolism (Vachon & del Giorgio, 2014;Vollenweider, 1975;Zwart et al., 2017). We choose to use precipitation minus evapotranspiration as our explanatory variable rather than streamflow because of the more direct connection to changes in climate.
To examine the effect of air temperature on C flux and productivity response variables, we regressed the average annual temperature against the residuals of the best fit linear model between respective C fluxes and the cumulative annual difference in precipitation and evapotranspiration. We use the statistical software package R for all model runs and statistical analyses.

Results and Discussion
The role of lakes in the global carbon cycle and their ability to support food webs via primary productivity are expected to change under future climate, but predictions of the magnitude and direction of response are extremely limited. Our analysis showed that the divergence among climate projections for the six GCMs resulted in a diversity of possible outcomes for regional C flux and productivity response to future climate. Cross-scale interactions between local hydrologic and physiochemical characteristics and regional drivers dictated the magnitude and direction of lake-specific carbon flux and productivity responses to simulated changes in climate. Interestingly, lake productivity and carbon cycling processes (e.g., CO 2 emission and C burial) did not respond in the same way to regional climate drivers because local characteristics mediated lake responses to climate change. These cross-scale interactions highlight the need for simulation models that explicitly capture these effects. Below we discuss the most important regional drivers and local characteristics for predicting lake C flux and productivity in response to future climate projections for our modeled region and opportunities for expanding to other regions.

Drivers of Regional Lake C Flux and Productivity Response to Future Climate
Variation in simulated regional water balance (annual streamflow or precipitation minus evapotranspiration) dictated whether lakes collectively were to increase or decrease in CO 2 emissions and C burial in our modeled region (Figures 1 and 2). Wetter simulations increased lake C emissions because increased streamflow also increased inputs of DIC (maximum increase in DIC loading of 29% in wettest scenario, supporting information Table S2), which, once transported to the lakes, was emitted to the atmosphere as CO 2 (Figure 2a). Additionally, increased heterotrophic respiration from increased organic carbon load (OC; maximum increase in OC loading of 28% in wettest scenario, supporting information Table S2) and increased water temperatures resulted in increased C emissions. The wetter simulations decreased WRT (maximum decrease in median lake WRT of 105 days in wettest scenario, supporting information Table S2), and thus, the processing time of C so that less C was removed within lakes relative to the amount exported from the landscape (Figure 2d). The decreased percent C removed and increased C emissions in response to wetter simulations reflect the hydrologically mediated trade-off between the magnitude of C fluxes and processing time, as wetter simulations increased C loading while simultaneously decreased WRT ). Total regional lake C burial was also higher during wetter simulations due to increased particulate C load from the terrestrial landscape and higher lake gross primary production (GPP) stimulated by increased nutrient load (Figure 2e; supporting information Table S2), which sedimented out of the water column and was subsequently buried (Figure 2b). Despite increased C burial in wetter simulations, the difference between emissions and burial tended to be higher (Figure 2c), indicating that lakes in regions with wetter futures will likely accelerate climate change impacts, while lakes in drier regions will likely mitigate them. This result is consistent with previous large-scale modeling studies (Cardille et al., 2009).
In addition to changes in the regional water balance, variation in simulated air temperature was an important regional driver for predicting regional lake metabolic response to climate (GPP and NEP) but was not important for predicting changes in regional lake C emissions or burial (Figures 1 and 2). After accounting for annual precipitation minus evapotranspiration, average annual simulated air temperature explained the residual variation in regional GPP (p < 0.001, R 2 = 0.74), NEP (p < 0.001, R 2 = 0.95), and percent C removed (p < 0.001, R 2 = 0.83; supporting information Figure S2); however, this was not the case for regional C emissions, C burial, or the difference between C emissions and burial (all p values >0.4 for air temperature explaining residual variation; supporting information Figure S2). The direct effects of air temperature on organismal metabolism are well known and are typically positive effects (Kraemer, Chandra, et al., 2017;Yvon-Durocher et al., 2010). Indeed, after accounting for annual precipitation minus evapotranspiration, average annual air temperature was positively related to volumetric rates of GPP in the epilimnion for our modeled region, an indicator of organismal metabolism (mol C m −3 ; p = 0.032, R 2 = 0.35, r = 0.59; supporting information Figure S3).
Changes in air temperature can also affect lake stratification strength (Winslow et al., 2017), ice duration (Sharma et al., 2019), and evaporation (Wang et al., 2018), which can indirectly affect whole-lake primary productivity through changes in suitable habitat for phytoplankton to grow (Thackeray et al., 2008). These indirect effects of air temperature on phytoplankton habitat and direct effects of temperature on metabolism make predicting lake primary production under climate warming scenarios challenging (Kraemer et al., Figure 2. Projected regional change in total (a) lake CO 2 emissions, (b) lake C burial, (c) difference between C emissions and burial, (d) percent C removed within the lakes, (e) total gross primary production (GPP), and (f) total net ecosystem production (NEP) all as a function of the difference between annual precipitation and evapotranspiration for historic (blue circle), 2050s (orange squares), and 2080s (red triangles). The numbers inside the shapes of the projected lake C dynamics correspond to the six GCMs, which are numbered by their performance for the modeled region, 1 being the best 6 the worst as calculated by Byun and Hamlet (2018); supporting information Table S1). The solid black lines represent the linear best fit line, and the gray dashed line represent no change in regional C flux or productivity when compared to historic model run.
2017; Kraemer, Chandra, et al., 2017). Indeed, after accounting for annual precipitation minus evapotranspiration, air temperature was negatively correlated to total regional GPP (mol C; p < 0.001, R 2 = 0.74, r = −0.87 for air temperature explaining residual variation in total GPP; supporting information Figure S2e) through the negative effects of air temperature on suitable habitat for phytoplankton production, even though there was a positive temperature effect on volumetric rates of GPP (supporting information Figure S3b). Elevated air temperature increased lake surface evaporation both by lengthening the open-water period and increasing evaporation rates, thus reducing lake water volume (Sharma et al., 2019;Wang et al., 2018;Figures 1 and S4c). Increased evaporative losses of water also concentrated light attenuating constituents such as DOC, which reduces heat penetration and shrinks the upper mixed layer where most phytoplankton production occurs (Houser, 2006). Our model results suggest that regions with drier and warmer futures will have decreased lake primary production through less watershed nutrient input and reduction in lake water volume, which could have major consequences for freshwater food webs, fisheries, and reliable protein sources for some regions (Gaeta et al., 2014;Lynch et al., 2016).
Regional changes in water balance and air temperature were also important for predicting changes in NEP (Figure 2f). Wetter simulations decreased NEP because of increased organic C flux from the surrounding landscape (supporting information Table S2), which increased heterotrophic activity more than the increased autotrophic activity (Figures 2e and 2f). This pattern has been demonstrated previously for temperate and boreal regions as increased runoff decreased lake NEP despite increased GPP for both individual storm events and inter-annual precipitation regimes (Ojala et al., 2011;Vachon & del Giorgio, 2014;Zwart et al., 2017). Heterotrophy is also more sensitive to changes in temperature than autotrophy (Yvon-Durocher et al., 2010), and our results show that simulated air temperature had a significant negative effect, more heterotrophy, on total NEP after accounting for annual precipitation minus evapotranspiration (Figures 2f and S2f; p < 0.001, R 2 = 0.96, r = −0.98 for air temperature explaining residual variation in total NEP).

Local Characteristics Mediate Lake-Specific Response to Future Climate
Local lake hydrologic characteristics, in particular the fraction of lake hydrologic export as evaporation (FHEE), mediated lake-specific C flux responses to simulated future regional water balance (Figure 3). For example, lakes with higher FHEE (lakes more hydrologically isolated from the landscape) tended to increase in C emissions in both wet and dry simulations, while change in emissions for lakes with lower FHEE (lakes more hydrologically connected to the landscape) was highly dependent on simulated regional water balance (Figures 3a and S4a). Eighty-nine percent of the most hydrologically isolated lakes (the highest 10% FHEE lakes for each simulation) increased in C emissions across all climate simulations, while only 42% of the most hydrologically connected lakes (lowest 10% FHEE) increased in C emissions ( Figure 3a). This contrast in change in CO 2 emissions was driven by the dominant process supporting CO 2 supersaturation for each lake type, as hydrologically connected lakes' CO 2 emissions were driven by externally loaded DIC (mean DIC load to DIC internally produced = 11.4), while hydrologically isolated lakes' CO 2 emissions were driven by heterotrophic activity (mean DIC load to DIC internally produced = 0.27). In drier climate simulations, externally loaded DIC was reduced due to decreased hydrologic inflows, which disproportionately decreased CO 2 emissions for the hydrologically connected lakes, while heterotrophic activity often increased due to elevated water temperatures across all climate simulations which increased lake CO 2 emissions, especially for the more hydrologically isolated lakes (supporting information Figure S4a). Differential support of lake CO 2 supersaturation via hydrologic pathways has been shown previously (Bogard & Giorgio, 2016;Maberly et al., 2013;McDonald et al., 2013;Stets et al., 2009;Vachon & del Giorgio, 2014;Wilkinson et al., 2016); thus, understanding lake hydrologic characteristics (e.g., FHEE), and potential future changes in these characteristics, can help inform how individual lakes, and lakes in aggregate, will respond to broad-scale changes in future climate.
Both lake hydrologic characteristics and historic lake DOC concentration mediated lake-specific volumetric GPP response to simulated climate (Figure 4). DOC is often considered a master variable for regulating lake ecosystem function (Solomon et al., 2015), and our model showed that changes in DOC are strongly related to changes in volumetric GPP as lake-specific historic DOC concentration (average from 1980 to 2010) dictated whether there was a positive or negative relationship between the change in GPP and change in DOC (Figure 4b). Lakes with low historic DOC concentration (<10 mg C L −1 ) generally exhibited a positive relationship between volumetric GPP and DOC concentration change, while lakes with high historic DOC concentration (>15 mg C L −1 ) exhibited a negative relationship between volumetric GPP and DOC concentration change. The latter effect occurs because DOC typically originates from allochthonous Figure 3. Projected lake-specific change in (a) lake CO 2 emissions, (b) lake C burial, (c) the difference between lake C emissions and C burial, (d) percent C removed within the lake, (e) total GPP, and (f) total NEP all as a function of the fraction of hydrologic export as evaporation (FHEE). Each point represents a lake in a given future climate scenario, and the color of the point represents the projected average difference in precipitation and evapotranspiration for the region. Since we used 12 future scenarios, each lake appears as a point 12 times in each panel. The 12 colored lines represent each scenario's LOESS best fit line, and the black dotted line represents no change from historic model run. The marginal plots represent the density of lakes for a given response variable (y axes) and predictor variable (FHEE, x axis) for each scenario, again colored by the scenario's difference in precipitation and evapotranspiration. sources (Wilkinson et al., 2013) and negatively impacts phytoplankton production through its light attenuating properties at high concentrations (Karlsson et al., 2009). Nutrients often accompany DOC and are coexported from the landscape to lakes, which positively impacts phytoplankton production (Corman et al., 2018). Thus, at low DOC concentrations, increasing DOC simultaneously increases light limitation but more rapidly reduces nutrient limitation for primary production, while at higher DOC concentrations, increasing DOC concentration reduces GPP as light limitation outweighs the release from nutrient limitation (Seekell et al., 2015). This differing response of lakes with low and high DOC loading has been documented previously both empirically (Seekell et al., 2015) and with simple productivity models Vasconcelos et al., 2018).
Changes in lake DOC concentration and consequently its effects on volumetric GPP were influenced by lake hydrologic characteristics and regional water balance (Figure 4a). Seventy-four percent of lakes increased in DOC across all future climate simulations, and the largest increases in DOC concentration occurred for the most hydrologically isolated lakes (higher FHEE) in drier simulations (Figure 4a). Browning of inland waters through increased DOC has been documented over much of the Northern Hemisphere (Monteith et al., 2007), and our model results show that this browning can be induced by changes in regional climate, especially through evapoconcentration in drier scenarios. Elevated DOC can simultaneously stimulate phytoplankton production through increased associated nutrients (Seekell et al., 2015) and shrink the upper mixed layer of lakes due to reduced heat penetration from increased light attenuation (Houser, 2006). Indeed, 84% of the lakes that increased in DOC concentration also increased in volumetric GPP (mol C m −3 ; Figure 4b); however, shrinking habitat area offset the increased volumetric GPP resulting in only 33% of elevated DOC lakes that increased in total GPP (mol C; supporting information Figure S4c).

Expanding to Other Regions and Integrated Modeling
Our model results showed that changes in lake C fluxes and productivity were strongly related to projected changes in regional water balance, and we suspect that this relationship would generally hold for other regions. However, precipitation and temperature projections are heterogenous both regionally and globally (Byun & Hamlet, 2018;Schewe et al., 2014), indicating that not all regions are likely to respond similarly to future climate in terms of lake C flux responses. For example, in summer, the projections for our modeled region show more drier scenarios than wetter scenarios when compared to historic climate ( Figure 1). Likewise, relatively small increases in streamflow are projected for basins in other northern latitude areas

10.1029/2019GL083478
Geophysical Research Letters of the Midwest and Great Lakes region that have had considerable snow cover historically, while basins in southern latitudes are more likely to experience substantial increases in streamflow (Byun et al., 2019). Therefore, we might expect southern areas of the Midwest and Great Lakes region to increase in lake C fluxes and GPP under future climate, while more northern areas might decrease in lake C fluxes and GPP since drier conditions are most probable. However, given the cross-scale interactions that dictated lake responses to climate in our analysis, integrated process-based models should be run at regional scales to better inform our opinion of how lakes in these different regions may respond to projected changes in climate.
Our analysis focused on the impact of projected climate change on lake C fluxes and primary productivity. There are likely future changes to terrestrial biogeochemistry and anthropogenic changes that will also alter regional aquatic C fluxes that we did not incorporate in our model. For example, land use change can alter inflowing constituent concentrations such as nutrients and particulate carbon, which could have impacts on lake biogeochemistry. For example, if portions of our modeled region were converted from forested land to agriculture, nutrient application would wash off the agricultural fields and into freshwater systems, stimulating algal growth (Carpenter et al., 1998). Additionally, lake CO 2 emissions have been shown to be positively influenced by terrestrial productivity (Maberly et al., 2013), indicating that interactions between future climate and terrestrial productivity may influence lake C fluxes. Although we focused on the effects of climate change on lake C fluxes and productivity via changes to the regional hydrologic cycle and lake physics, broad-scale integrated modeling that incorporates human decision making and terrestrial and aquatic interactions is an important and challenging avenue for future modeling research (Carpenter et al., 2015;Gil et al., 2018;Rose et al., 2017).

Conclusions
Lakes contribute significantly to global biogeochemical cycling, provide habitat for fish and wildlife, and offer a wide range of ecosystem services to humans. Anticipating how these important ecosystems may change in future climates is imperative. Despite high regional and lake-specific variability in lake C flux and productivity (Figures 1 and 2), we show that regional lake C flux and productivity were strongly related to regional climate drivers including water balance and air temperature and lake-specific changes to these drivers were mediated by local characteristics including lake hydrologic characteristics and historic DOC concentration (Figures 2-4). Generally, in regions that become wetter, lakes will likely increase in net C export (increased C emissions minus C burial) and the change in lake primary productivity will depend on the magnitude of air temperature increase. In regions that become drier, lakes will likely decrease in net C export (decreased C emissions minus C burial); however, lake primary productivity will be greatly reduced with consequences for the rest of the aquatic food web. Given the heterogenous air temperature and precipitation changes both globally and regionally (Byun & Hamlet, 2018;Schewe et al., 2014), future lake contributions to global C cycles and primary productivity remain uncertain; however, process-based models that account for these cross-scale interacting effects of regional drivers and local characteristics on lake C fluxes and productivity provide insight into how these complex systems may respond to projected changes in climate.