Reviewer 1

Major Changes To Methodology (Sensitivity, Calibration, and Uncertainty)

Based on reviewer comments, we added a sensitivity analysis of the model parameters. We used the global sensitivity method from  Morris (1991) to investigate the sensitivities of model output variables on each model parameter. The sensitivity analysis showed that there were nine parameters to which the model was consistently sensitive across the six study lakes. This group included the ratio of DOC and POC produced from NPP (C_npp), the maximum daily productivity parameter (Pmax), the inflow concentration of recalcitrant POC (POCR_inflow), the setting velocity of recalcitrant POC (K_POCR), the temperature fitting coefficients for productivity and respiration (θ_npp, θ_resp), the slope of the irradiance/productivity curve (IP), the sediment respiration flux (Rsed), and the respiration rate of DOCL (Resp_DOCL). 

We chose a subset of the nine parameters to include in the uncertainty analysis based on the following justifications. The model results showed that recalcitrant substrates are of lesser importance for lake metabolism dynamics, so we chose not to further investigate the uncertainty of the POCR_inflow and K_POCR parameters. The Pmax and IP parameters are directly correlated, so we chose to remove Pmax from further uncertainty considerations. The θ_npp and θ_resp parameters act as substitutes for water temperature, a well-known “master variable” in water quality modeling, and directly reflect seasonality in the model. Therefore,we chose to omit these parameters for further uncertainty calculations. The final subset of parameters for uncertainty analysis consisted of C_npp, Resp_DOCL, R_sed, and IP. Of the four parameters, we felt C_npp was best constrained by the literature. To reduce the number of parameters estimated in the calibration process we restricted the automated constrained parameter search to the remaining three.

To clarify the text regarding model calibration we rewrote that section of the methods to better integrate the sensitivity analysis. There are three types of parameters – constants, manually calibrated based on literature, and constrained search based on the previous sensitivity analysis and literature ranges. Parameters with low sensitivity on model outcomes were manually calibrated. On the contrary, to account for the high sensitivity of IP, Rsed, and Resp_DOCL, we ran an automatic constrained search across the joint parameter space to identify values that led to the best model fit.

Sensitivity guided uncertainty analysis. To quantify uncertainty around model predictions, we sampled IP, Rsed, and Resp_DOCL simultaneously from uniform distributions defined by literature values (Table 3). We plotted the 2.5 and 97.5% quantiles for these distributions and included them in the time series plots (Fig 2, 3, 4, SI Fig 1,2,3).

Specific Comments: 

Line 79: It would be helpful to provide more detail on “external forcings” here. Which variables are you most considering? Hydrology, climate, nutrient inputs?

Reply by authors: We have updated the sentence to explain this in more detail: 

“An increase in the prevalence of hypolimnetic anoxia and associated water quality degradation in temperate lakes indicates the need to better understand how lake ecological processes interact with external forcings, such as hydrology and nutrient inputs, to lead towards the development of anoxia (Jane et al. 2021)”

Line 88-96: Allochthonous OC can also impact DO dynamics through increasing thermal stratification due to its high chromophoricity. Is this process included in the model? Can you disentangle allochthonous OC influences from the effect of increased stratification versus directly fueling microbial respiration?

Reply by authors: Stratification dynamics in the model are estimated from previously modeled temperature profiles for the study lakes. In the previously modeled temperature profiles, water characteristics controlling thermal stratification, such as chromophoric organic matter, would have been subsumed in other parameters during model fitting. The time dynamics of chromophoricity of allochthonous OC is not considered for thermal stratification in our model. However, the light extinction of DOC (LEC_DOC) is a parameter that we use in the model to calculate the total light extinction coefficient (K_LEC), which in turn is used to control NPP. We assume that the majority of DOC in a lake is from allochthonous sources, so for lakes with higher allochthonous loads we increase the LEC_DOC parameter to represent the effect that darker (newer) DOC has on water column light and the NPP of a lake (SI Table 5).

Figure 1: To improve clarity and simplicity of the model conceptual figure, consider removing the dashed boxes of TP and IC, as they are not state variables tracked in the model. Specifically, since IC is considered non-limiting and not a driver in any equation, the visualization would benefit from its removal. TP could be included by representing primary productivity as a function dependent on TP, ie. primary productivity(TP).

Reply by authors: We agree and have made these changes.

Table 3, eq. 20: Do recalcitrant and labile OC have different chromophoricity? They are treated similarly in their effect on light extinction coefficient in the model. Did you consider weighting recalcitrant OC more due to high aromaticity or light attenuating compounds? Would that change model dynamics?

Reply by authors: We do not account for chromophoricity differences between recalcitrant and labile OC. Given its labile nature, autochthonous DOC rarely accounts for more than 10-20% of the total DOC pool. DOC observed in a lake is mostly from allochthonous sources that are more recalcitrant in nature. While we can imagine very interesting short-term dynamics in water clarity as a function of, e.g., storm events, our focus in this manuscript was on longer term dynamics. 

What about all the other parameters that were manually fit in Table 4. How were those fit and was any sensitivity analysis conducted?

Reply by authors: Please see our top reply regarding Major Changes to Methodology

Line 394: Free parameters were fit manually across their ranges, and it seems parameter values chosen were often at the extremes of the ranges. They also are strongly segregated across the northern and southern lakes, which may be driving many of the differences. How much did parameter value influence model fit? A sensitivity analysis here is needed, or atleast reporting the model residuals with different parameter values. Additionally, did you consider extending the ranges since chosen values were often at the extremes? 

Reply by authors: Please see our top reply regarding Major Changes to Methodology

Lines 473-476: The authors report differences in the contributions of OC fates in the overall budget across the study lakes. Are these statistically significant differences? It would help interpretation to determine this, as the error bars on Figure 5 seem to be overlapping across lakes, but it is hard to compare across lake panels.

Reply by authors: We compared the contributions of OC sources and fates across the study lakes using the Wilcoxon statistical test for significance. Please see the attached table (SI Table 6) for these results. Note, any p-value less than 0.05 was considered significant. All non-significant differences between lakes are indicated by “NS”. We found similarities in the flux values between ME and MO as well as BM and SP. This helps highlight that OC fluxes for lakes within the same trophic status were more similar and that OC fluxes for lakes 

across the trophic gradient were less similar.

Lines 580-583: What about anaerobic carbon metabolism? How would anaerobic metabolism influence sediment OC pools? The surplus of auto OC could be a direct reflection of sediment respiration nearing zero when oxygen is depleted. Previous research shows relationships between anaerobic OC degradation and autochthony, such that the OC could be mineralized as CH4, and there is a growing body of work on the importance of anaerobic metabolism for overall carbon cycling and OC burial in stratified lakes. I think it is worth noting that this is not considered and considering how it might change model predictions.

Reply by authors: We neglected anaerobic metabolism in the current modeling work to focus intensively on the relationship between carbon cycling and oxygen depletion. Nonetheless, as the reviewer has pointed out, anaerobic mineralization of organic carbon is an important biogeochemical process. Under low DO concentrations, anaerobic mineralization can account for up to > 92 % of organic carbon respiration mostly through methanogenesis (Maerki et al. 2009). Further methane oxidation can be an important sink for oxygen. We added a clarification to the Methods section 2.3.4 (Internal Lake Metabolism Fluxes):

“Please note that we did not include anaerobic carbon metabolism in our modeling approach and discuss potential shortcomings in the discussion section”

Further, we added this paragraph to the discussion:

“Anaerobic mineralization of organic carbon is an important biogeochemical process that can be the main carbon sink through methanogenesis (Maerki et al. 2009) and which was not incorporated into our modeling scheme. We envision future metabolism studies to focus on incorporating these processes into their modeling schemes, which will potentially decrease the amount of OC buried annually. As the main source for sediment OC in our model is autochthonous carbon, we envision that anaerobic mineralization would (a) decrease autochthonous OC in the hypolimnion, and (b) decrease the overall OC sediment pool. Nonetheless, loss of OC from the system through anaerobic respiration would primarily happen through methane  production. Hart (2017) highlighted that there is little evidence of a diffusive flux of methane out of Lake Mendota, probably due to its oxidation in the water column. Methane oxidation would be subsumed under aerobic respiration in our modeling approach. Only methane ebullition is not accounted for in our model. Although we believe that ebullition is not a substantial portion of the lake’s carbon mass budget, that remains to be studied more carefully.  As the model accounts for DO consumption through calibration, the overall flux would not change even if we link DO consumption to methane oxidation, only the process description would be more realistic.“

Technical Corrections:

Line 176-192: Where was hydrology data from Lake MO obtained from? This lake is missing from the description.

Reply by authors: A description for Lake MO hydrology can be found in the manuscript (L264-266): 

“ME is the predominant hydrologic source for MO (Lathrop & Carpenter 2014), thus, MO inflow is assumed to be equal to ME outflow, and MO outflow is assumed to be equal to MO inflow.”

Table 3, Eq. 1: Missing open parentheses in front of NPP.

Reply by authors: Thank you for catching this. We have fixed this and it will be updated for the new version of the manuscript.

Additional References:

Maerki, Martin, Müller, Beat, Dinkel, Christian, Wehrli, Bernhard, (2009), Mineralization pathways in lake sediments with different oxygen and organic carbon supply, Limnology and Oceanography, 54, doi: 10.4319/lo.2009.54.2.0428.

Morris, M. D.: Factorial Sampling Plans for Preliminary Computational Experiments, Technometrics, 33, 161–174, https://doi.org/10.1080/00401706.1991.10484804, 1991.

Reviewer 2

Major Changes To Methodology (Sensitivity, Calibration, and Uncertainty)

Based on reviewer comments, we added a sensitivity analysis of the model parameters. We used the global sensitivity method from  Morris (1991) to investigate the sensitivities of model output variables on each model parameter. The sensitivity analysis showed that there were nine parameters to which the model was consistently sensitive across the six study lakes. This group included the ratio of DOC and POC produced from NPP (C_npp), the maximum daily productivity parameter (Pmax), the inflow concentration of recalcitrant POC (POCR_inflow), the setting velocity of recalcitrant POC (K_POCR), the temperature fitting coefficients for productivity and respiration (θ_npp, θ_resp), the slope of the irradiance/productivity curve (IP), the sediment respiration flux (Rsed), and the respiration rate of DOCL (Resp_DOCL). 

We chose a subset of the nine parameters to include in the uncertainty analysis based on the following justifications. The model results showed that recalcitrant substrates are of lesser importance for lake metabolism dynamics, so we chose not to further investigate the uncertainty of the POCR_inflow and K_POCR parameters. The Pmax and IP parameters are directly correlated, so we chose to remove Pmax from further uncertainty considerations. The θ_npp and θ_resp parameters act as substitutes for water temperature, a well-known “master variable” in water quality modeling, and directly reflect seasonality in the model. Therefore,we chose to omit these parameters for further uncertainty calculations. The final subset of parameters for uncertainty analysis consisted of C_npp, Resp_DOCL, R_sed, and IP. Of the four parameters, we felt C_npp was best constrained by the literature. To reduce the number of parameters estimated in the calibration process we restricted the automated constrained parameter search to the remaining three.

To clarify the text regarding model calibration we rewrote that section of the methods to better integrate the sensitivity analysis. There are three types of parameters – constants, manually calibrated based on literature, and constrained search based on the previous sensitivity analysis and literature ranges. Parameters with low sensitivity on model outcomes were manually calibrated. On the contrary, to account for the high sensitivity of IP, Rsed, and Resp_DOCL, we ran an automatic constrained search across the joint parameter space to identify values that led to the best model fit.

Sensitivity guided uncertainty analysis. To quantify uncertainty around model predictions, we sampled IP, Rsed, and Resp_DOCL simultaneously from uniform distributions defined by literature values (Table 3). We plotted the 2.5 and 97.5% quantiles for these distributions and included them in the time series plots (Fig. 2, 3, 4, SI Fig. 1,2,3).

General Comments:

The introduction struggles to highlight the need to connect the C and O cycles in lakes. Many sentences are vague and poorly informative, leaving the reader to check by themselves what is really said in the cited references. For instance, at different places, trophic status, hydrology, and morphometry are cited as all acting/influencing deep hypoxia (l74,l92,l100, l107, l109). But how they are mechanistically acting upon hypoxia remains very vague, and readers should go and check the references. Instead of mentioning that factors 'influence" hypoxia, I would find it much more informative to explain how they do, providing more informative facts. I would suggest rewriting and refocusing the introduction, enriching the conceptual background. The main achievement of this paper is to show that heterotrophy does not necessarily generate greater hypoxia. There is a body of research relating terrestrial organic matter inputs to lake respiration; then, an easy set of hypotheses would be to drive from this. With a greater allochthony and heterotrophy, greater respiration could shift the oxygen balance towards depletion. Then, I suggest focusing on the processes acting upon hypolimnetic oxygen, namely, those by which oxygen is consumed and then relates to carbon sources, and those by which oxygen is renewed. It would be finally easier, from processes, to introduce how lake hydrology and morphometry can affect both consumption and renewal and finally modulate what to expect regarding oxygen depletion.

Reply by authors: The reviewer makes excellent suggestions for reformulating the Introduction to the manuscript. In keeping with the suggestions, we shift the emphasis to the possible role of autochthony to hypolimnetic respiration as an advancement to the oft-cited paradigm of net heterotrophy being the norm for lakes, with the implication that gradients of allochthony should explain gradients of hypolimnetic oxygen depletion. We provide the following topic sentences for a new Introduction, and we hope to have the opportunity to share the entire Introduction should we move to the next stage of the publication process.

1. Hypolimnetic oxygen depletion is a persistent and global phenomenon that degrades lake ecosystems services.
2. Allochthonous organic carbon (OC) loading to lakes explains the prevalence of negative net ecosystem production (i.e., net heterotrophy) and provides substrate for hypolimnetic oxygen depletion. 
3. The contributions of OC from autochthony to hypolimnetic oxygen depletion may be important as well, despite its low concentrations relative to that of allochthonous OC in many lakes.
4. Understanding the relative importance of autochthony and allochthony to hypolimnetic oxygen depletion requires consideration of a number of physical and biological processes controlling oxygen sources and sinks in lakes.
5. The availability of long-term observational data combined with process-based models provides an opportunity to investigate OC sources and their control over the dynamics of lake DO across multiple time scales.
6. In this study, our goal is to investigate OC source contributions to lake carbon cycling and hypolimnetic oxygen depletion. We address the following questions: 
1. What are the dominant sources of organic carbon that contribute to hypolimnetic oxygen depletion, and how do their contributions differ across a group of diverse lakes over two decades? 
2. How do lake trophic state, hydrology, and morphometry influence the processing and fates of organic carbon loads in ways that affect hypolimnetic dissolved oxygen?

I still have not understood how allochthonous loadings were computed. Discharges come essentially from a hydrological model, but what about concentrations?

Reply by authors: Our process for computing allochthonous loads for the study lakes is described in L255-259 and L267-270. However, we believe that the reader would benefit from a more detailed description and have updated Methods section 2.3.2 (External Lake and Environmental Physical Fluxes) with the following information:

“ME allochthonous loads are based on observed allochthonous DOC and POC concentration ranges found in Hart et al. (2017) and were verified against the modeled OC equilibrium for the lake. MO inflow concentrations for DOC and POC were taken from the surface water concentrations of ME and were verified against load estimates in McCullough et al., 2018. Northern lake total OC inflow concentrations were taken from estimated annual loads found in Hanson et al., 2014. Specifically, DOC inflow concentrations were calibrated from winter in-lake observational data. POC inflow concentrations were calibrated from the remaining OC inflow budget for lakes.”

There is also sometimes confusion between what is a flux and what is a rate. For instance, l303, "sediment respiration for the hypolimnion [...] is a constant daily rate", rate should be included after respiration.but l 332, "the respiration rates are also scaled by oxygen availability", here rate should be replaced by "flux"

Response by authors: This was an oversight on our part, and thank you for bringing it to our attention. We have fixed this in the manuscript.

l390, acronyms for r_sed and R_DOCL are different (Resp_DOCL and Resp_sed)

Response by authors: We have changed these acronyms to make them consistent throughout the manuscript.

Specific Comments:

R: l69-74. It is surprising to start the introduction by mentioning that hypoxia can lead to increased OC release from the sediment. At the same time, this process is not included in the very model aimed at linking the O and C cycle. The study of drivers of hypoxia in lakes could be better motivated.

Response by authors: We have removed this sentence from the introduction. 

Table 3: r_rate not defined

Response by authors: We have added the following text to the Table 2 description: 

“The term (r_rate) is included in Eq. 13 to represent the respiration rates of the different OC pools. It is included to simplify the table of equationsTerms not defined here are included in Table 3.”

Results: RMSE is the only metric used to assess the goodness of fit. However, it poorly captures the model's ability to reproduce the seasonal dynamics, while this is the critical requirement for the model (as it is reinitialized each Jan 1st). Does it make sense to add another metric?

Response by authors: In addition to RMSE, we have calculated the Nash-Sutcliffe model efficiency coefficient (NSE) and the Kling-Gupta Efficiency (KGE) for the model state variables. These values are presented in SI Table 7.

L615- "It has been shown that POC respiration contributes substantially to hypolimnetic DO depletion (Jenny et al. 2016)," This quote needs to be checked. I do not recall Jenny et al quantifying the respiration of POC.

Response by authors: The reviewer has identified an oversight in our referencing the literature. We have updated the quote:

“POC respiration can contribute substantially to hypolimnetic DO depletion in both lakes and reservoirs (Beutel, 2003)”

L658-662. Is the long-term effect of allochthony testable from a model that uses annual reinitialization?

Response by authors: The model states are initialized at the first time step of the model run, and are not reinitialized at any point. We do include annual stratification and mixing dynamics, which are estimated using temperature profiles to calculate water column buoyancy. We have updated section 2.3 of the manuscript (“The Model”) to clarify this point. See below:

“Throughout each year, the model tracks state variables and fluxes in the lake for each day (Fig. 1). These state variables are initialized at the first time step of the model and include DO and the labile and recalcitrant components of particulate organic carbon (POC) and dissolved organic carbon (DOC).

Additional References:

Morris, M. D.: Factorial Sampling Plans for Preliminary Computational Experiments, Technometrics, 33, 161–174, https://doi.org/10.1080/00401706.1991.10484804, 1991.

Beutel, Marc (2003) Hypolimnetic Anoxia and Sediment Oxygen Demand in California Drinking Water Reservoirs, Lake and Reservoir Management, 19:3, 208-221, DOI: 10.1080/07438140309354086

Reviewer 3

Major Changes To Methodology (Sensitivity, Calibration, and Uncertainty)

Based on reviewer comments, we added a sensitivity analysis of the model parameters. We used the global sensitivity method from  Morris (1991) to investigate the sensitivities of model output variables on each model parameter. The sensitivity analysis showed that there were nine parameters to which the model was consistently sensitive across the six study lakes. This group included the ratio of DOC and POC produced from NPP (C_npp), the maximum daily productivity parameter (Pmax), the inflow concentration of recalcitrant POC (POCR_inflow), the setting velocity of recalcitrant POC (K_POCR), the temperature fitting coefficients for productivity and respiration (θ_npp, θ_resp), the slope of the irradiance/productivity curve (IP), the sediment respiration flux (Rsed), and the respiration rate of DOCL (Resp_DOCL). 

We chose a subset of the nine parameters to include in the uncertainty analysis based on the following justifications. The model results showed that recalcitrant substrates are of lesser importance for lake metabolism dynamics, so we chose not to further investigate the uncertainty of the POCR_inflow and K_POCR parameters. The Pmax and IP parameters are directly correlated, so we chose to remove Pmax from further uncertainty considerations. The θ_npp and θ_resp parameters act as substitutes for water temperature, a well-known “master variable” in water quality modeling, and directly reflect seasonality in the model. Therefore,we chose to omit these parameters for further uncertainty calculations. The final subset of parameters for uncertainty analysis consisted of C_npp, Resp_DOCL, R_sed, and IP. Of the four parameters, we felt C_npp was best constrained by the literature. To reduce the number of parameters estimated in the calibration process we restricted the automated constrained parameter search to the remaining three.

To clarify the text regarding model calibration we rewrote that section of the methods to better integrate the sensitivity analysis. There are three types of parameters – constants, manually calibrated based on literature, and constrained search based on the previous sensitivity analysis and literature ranges. Parameters with low sensitivity on model outcomes were manually calibrated. On the contrary, to account for the high sensitivity of IP, Rsed, and Resp_DOCL, we ran an automatic constrained search across the joint parameter space to identify values that led to the best model fit.

Sensitivity guided uncertainty analysis. To quantify uncertainty around model predictions, we sampled IP, Rsed, and Resp_DOCL simultaneously from uniform distributions defined by literature values (Table 3). We plotted the 2.5 and 97.5% quantiles for these distributions and included them in the time series plots (Fig. 2, 3, 4, SI Fig. 1,2,3).

General Comments: 

I could not determine whether the metabolism model was specifically developed for this analysis or whether it was based on some previously published model, as no citations were offered. If the model is specific to this analysis, citations should be provided that justify different modeling choices. 

Reply by authors: The model was coded specifically for this analysis; however, many of the assumptions around model complexity and mathematical formulations were borrowed from literature cited. We chose to develop our own process-based model rather than use an existing model (such as GLM, Simstrat, etc.) so that we could simulate and measure the specific metabolism fluxes related to our study questions. As the reviewer points out, there are some important differences. We clarify these points through additional annotations and citations in the manuscript, as well as an additional paragraph written to clarify the issue. 

I couldn’t find a citation for the model linking TP to NPP (Equation 12). Where did this come from?

Reply by authors: Total phosphorus concentration in a layer is from observational data for each lake interpolated to the daily time scale. The interpolated values are then normalized for each individual lake to drive NPP. These values are standard-normal transformed for the entire time series; thus, the mean and variance of phosphorus is subsumed in the estimates of the IP and Pmax parameters. The time dynamics of normalized phosphorus concentrations are retained to represent seasonal P dynamics in the lake.

What is Pmax?

Reply by authors: The Pmax parameter is the maximum daily productivity for a lake and is based on trophic specific ranges provided in Wetzel (2001). These values are provided in units of [mgC/m2/day], and subsume lake-specific nutrient concentrations. We convert these values to units of [gC/m3/day] for use in our model.

How was this DOC allocated to labile and recalcitrant fractions? I’m guessing that it was assigned entirely to the recalcitrant fraction, but no information is provided.

Reply by authors: Allochthonus DOC loads were assigned entirely as recalcitrant substrates. To clarify this point for the reader, we have included the following description in Methods section 2.3.2 (External Lake and Environmental Physical Fluxes): 

“For the northern lakes (TR, AL, BM, SP), we assume that allochthonous OC loads consist of entirely recalcitrant substrates and calibrate these and our recalcitrant OC export values from Hanson et al. (2014).”

Also, the presence of a INDOCL term in Equation (2) suggests that at least some allochthonous load is designated as labile.

Reply by authors: We assume that allochthonous DOC loading is a recalcitrant substrate, with the exception of Lake Monona (MO) where external OC loads are estimated from Lake Mendota (ME). Given this relationship, we account for labile POC and DOC loading for MO and include additional terms for INDOCL and INPOCL for Equation 2 and 4 in Table 2, respectively.

Recalcitrant vs. labile OC. The focus of this paper is to determine the relative contributions of allochthonous and autochthonous OC to oxygen depletion, and the lability of OC is a crucial determinant of these contributions. As stated by the authors, autochthonous OC tends to be more labile than allochthonous OC, but as far as I can tell, all allochthonous OC is assigned to the recalcitrant fraction and all autochthonous OC is assigned to the labile fraction. If this is so, then isn’t the conclusion that most oxygen depletion is driven by autochthonous OC self-evident? Why bother running the model?

Reply by authors: For recalcitrant vs. labile OC above: The total respiratory flux is the product of the respective labile and recalcitrant pool sizes and their decay rates, the available oxygen, and the ambient temperature. The oxygen demand, therefore, depends very much on the OC pool sizes, which vary by seasons and differ between thermal strata. Based on ecosystem observations alone, the pool size of the recalcitrant fraction of OC is higher in many lakes, suggesting that ecosystem respiration scales with allochthony. Without quantifying primary production and its fate, we have no way to account for its contribution to total respiration. How OC moves through the lake matters, as well. For example, it is conceivable that primary production in the epilimnion could have been offset by epilimnetic respiration, thus eliminating the autochthonous pool before it had a chance to contribute to hypolimnetic respiration. However, through modeling, we were able to quantify the autochthonous OC production and its export from the epilimnion to the hypolimnion and its total contribution to the organic carbon cycle, including the effects on oxygen consumption and the seasonality of those dynamics. These ideas are emphasized in the rewriting of the Introduction and in additional points made in the Discussion.

Model uncertainty. More work is needed to quantify the uncertainty in the model. The uncertainty in the key components of the model (estimating NPP and external OC loads) is likely large, but no information is offered as to how this uncertainty would affect model results. Assuming that the uncertainty of all model parameters is 20% of the mean value is also a gross simplification.

Reply by authors: Please see our top reply regarding Major Changes to Methodology.

Line-Specific Comments:

Table 1: The authors state that the sampled lakes cover a broad range of allochthonous loads, but the range of DOC concentrations is narrow relative to concentrations observed at the continental scale (see for example, values reported in EPA’s National Lakes Assessment. It would be good to place the DOC concentrations in a broader context.

Reply by authors: The DOC concentrations provided in Table 1 cover a broad range of allochthonous loads associated with landscapes across the North Temperate Lakes district.

Line 190-192: The citation here points to documentation for the data but how exactly were nutrient concentrations used to computed PP?

Reply by authors: We compared our coefficients for primary production with those found in (Wetzel, 2001), which are in the following units [mgC/m2/day]. Note that such coefficients are not per unit phosphorus, but rather, subsume lake-specific nutrient concentrations. In order to do our comparison, we subsumed nutrient concentrations into our calculations of primary production. Our approach was to remove the mean of observed P and normalize by the P variance. Those two statistical features of P become subsumed in the estimates of IP and Pmax, but the time dynamics of normalized P are retained to represent seasonal P dynamics in the lake.

Line 204: The Hipsey et al. 2019 paper is cited repeatedly but in my reading, this paper only covers the hydrodynamic part of the lake model and says nothing about metabolism models. Is the citation incorrect, or is there a different paper that should be cited?

Reply by authors:The reviewer has identified an oversight in our referencing the literature. We have updated the source to Hipsey et al. 2022. 

Hipsey, M.R., (ed.). 2022. Modelling Aquatic Eco-Dynamics: Overview of the AED modular simulation platform. Zenodo. https://doi.org/10.5281/zenodo.6516222.

Line 257: What does the sentence here mean? It seems to state the same information as the previous sentence. (included below)

L257 – “We specifically use the allochthonous OC load values in this study to assist in the manual calibration of inflow recalcitrant POC and DOC concentrations for each lake”

Reply by authors: This sentence will be removed as it is redundant. 

Line 266 - 267: How is the calibration performed? The citation again just documents the source of the data. 

Reply by authors: This was an oversight on our part. Allochthonous OC concentrations were not calibrated. These values were taken from observed data in Hart et al. (2017) and we verified these values by back-calculating inflow concentrations based on the modeled OC equilibrium for the lake.

Line 319: So, all OC from NPP is assumed to be labile, right?

Reply by authors: Yes, this assumption is accurate. We have updated this sentence to reflect our assumptions.

“All OC derived from NPP is assumed to be labile and is split between particulate and dissolved OC production, with eighty percent produced as POC and twenty percent produced as DOC. This ratio was determined through model fitting and is similar to previously reported values (Hipsey et al. 2022).

Line 350: Is inorganic suspended sediment negligible in these lakes? Most Secchi computations would include a contribution from inorganic sediment.

Reply by authors: Due to the scarcity of long-term TSS data for these lakes, we have taken this into our model assumptions. We assume that TSS inorganic particulates are low in concentration and are a small percentage of the background water clarity. Given this, we do not include the impacts of inorganic TSS on Secchi in our model. 

Line 406: This approach for Secchi seems very ad hoc. Why are Secchi’s in northern lakes so variable? Higher sampling error? More temporal variability?

Reply by authors: We do not know why observed Secchi data are highly stochastic. We confirmed that characteristic in the observational data using time series decomposition. While there is an annual signal of about +/- 0.5 m, the noise component can be up to four times greater. As an example we have attached a time series decomposition for TR to this reply (SI Fig 7).

Line 442: The differences between model predictions of DOC and observed values are so great that it’s hard to argue that the model is accurately reproducing temporal changes in DOC. Only predictions TR and SP are close to observations and that seems more due to the fact that the magnitude of DOC temporal variations is small.

Reply by authors: We recognize that there are small biases in the long term values of the prediction of the doc. It is also the case that the variance is not fully reproduced by the model. We document the goodness of fit more fully in SI Table 7 (attached) for the benefit of the reader. 

Line 561: Good point. I’m struggling to determine whether the model described in this manuscript is accurate enough to support this statement.  A more realistic sensitivity and uncertainty analysis would help, particularly if that uncertainty was carried through to examine the effects on the budgets depicted in Figure 5.

Reply by authors: We do not assert that our model outcomes prove this point but rather support it. We did add a sensitivity analysis which is described in the top reply. Our analysis of the budgets helps support this (do statistical significance of OC sources as well and attach to this).

Additional References: 

Hipsey, M.R., (ed.). 2022. Modelling Aquatic Eco-Dynamics: Overview of the AED modular simulation platform. Zenodo. https://doi.org/10.5281/zenodo.6516222.

Morris, M. D.: Factorial Sampling Plans for Preliminary Computational Experiments, Technometrics, 33, 161–174, https://doi.org/10.1080/00401706.1991.10484804, 1991.