The manuscript addresses a significant and timely topic in cold-region hydrology. The methodology is sound, and the results provide a clear path toward improving operational streamflow forecasting. However, the manuscript requires a more critical discussion regarding the physical basis of the macropore parameterization and a deeper analysis of the spatial variability in model performance. Addressing the major comments regarding the physical justification and the sensitivity of the alpha_MP parameter will significantly enhance the impact and scientific rigor of the work. I look forward to seeing a revised version of this study.

 Major Comments：

1. The authors propose a non-explicit, threshold-based representation of macropores. While this is computationally efficient for operational forecasting, the manuscript lacks a robust discussion on the physical limitations of this simplification. Specifically, how does this approach account for the spatial heterogeneity of macropore networks at the grid-cell scale? A more thorough comparison with existing dual-domain or dual-permeability models (e.g., those cited in the introduction) is needed to justify why this simplified approach is sufficient for large-scale hydrological modeling.

2.  The calibration of alpha_MP is a central component of the study. However, the manuscript does not sufficiently address the potential for equifinality or the transferability of this parameter across different hydro-climatic regions. Are the optimal values of alpha_MP consistent across the diverse landscapes of the Great Lakes and Saint-Lawrence domain? A more detailed analysis of the spatial variability of the optimal alpha_MP and its relationship with soil/land-cover characteristics would significantly strengthen the paper.

3.  The results indicate that the Fr-MP configuration performs differently in agricultural areas compared to forested catchments. The authors attribute this to the interaction between tile drainage and frozen soil. I suggest a more in-depth analysis of this interaction. Does the current model structure adequately decouple the effects of macropores from those of anthropogenic drainage systems? This is crucial for the model's reliability in human-altered landscapes.

4.  While the authors evaluate the impact on TT and TD, the improvements are described as "slight" or "neutral." Given that the primary goal is to improve streamflow, the authors should clarify whether these meteorological improvements are statistically significant and if they have any meaningful impact on the broader numerical weather prediction (NWP) performance at ECCC.

 Minor Comments：

1.  Page 2, Line 40: The term "discrepancy" is used well here; ensure this is consistently linked to the scale-dependency of the processes throughout the discussion.

2.  Page 7, Eq. 9: Please clarify if alpha_MP is a static parameter or if it could be dynamic based on seasonal vegetation growth.

3.  Page 8, Figure 1: The conceptual diagram is excellent. Please ensure the legend clearly distinguishes between the "matrix flow" and "macropore flow" pathways in panel (c).

4.  Page 10, Section 2.4.1: The exclusion criteria for stations (NSE < 0 or KGE < -0.41) are reasonable, but please provide a brief justification for why these specific thresholds were chosen over others.

5.  Page 13, Line 355: The mention of the "perched aquifer" in Michigan is interesting. Could the authors provide a bit more detail on how this specific geological feature might be better represented in future iterations of the model?

6.  Page 16, Line 411: The observation that only 4 stations show a Delta KGE < -0.5 is a strong point. Please emphasize this in the abstract as a measure of the model's robustness.

7.  Page 23, Line 495: The authors mention that wetlands and meanders are not fully represented in Watroute. It would be beneficial to include a brief sentence on how the transition to the "Raven" routing scheme might specifically address these limitations.

8.  Clarity of Figures: In Figure 4 and Figure 5, the color scales are effective, but please ensure that the symbols for different drainage areas are easily distinguishable in the printed version.

9.  References: Ensure that all cited preprints (e.g., Bauer et al., 2025) are updated to their final publication status if available by the time of revision.

General comment:



This study presents a computationally efficient approach to improve infiltration under frozen soil conditions in a land surface model. The manuscript clearly demonstrates that the current soil freezing formulation leads to unrealistic runoff responses and that increasing infiltration can substantially improve streamflow simulations at large spatial scales. The topic is highly relevant, particularly in the context of operational hydrological forecasting, and the study provides a valuable and pragmatic contribution in this regard.

However, the proposed “macropore representation” is, in its current form, a conceptual parameterization that mimics the effect of preferential flow rather than explicitly representing macropore processes. The implementation effectively removes the ice-induced hydraulic impedance once a soil moisture threshold is exceeded and additionally modifies the runoff generation scheme. As such, the approach represents a combined structural modification of flow partitioning rather than a direct representation of macropore flow.

In its current form, the manuscript does not sufficiently disentangle these effects. A clearer separation of the impacts of the increased vertical hydraulic conductivity and the modified runoff generation would be necessary to properly attribute the improvements in model performance.

Furthermore, the physical interpretation of the results could be strengthened. In particular, the manuscript would benefit from a more explicit discussion of:

• the conceptual nature and limitations of the parameterization,
• the implications of representing subgrid-scale preferential flow using a bulk threshold formulation,
• the spatial variability in model performance, especially in relation to land use, soil texture, and the exclusion of stations from the evaluation.

Finally, the transferability of the calibrated parameter α_MP within and beyond the study domain is not addressed and should be discussed.

Overall, I consider the manuscript suitable for publication after revision. The proposed approach is particularly interesting from an operational perspective, as it provides a simple and efficient way to improve infiltration under frozen conditions. I am especially interested to see how this concept will evolve in future work toward more physically based representations of preferential flow processes.

Minor comments and technical notes:



L. 44: The sentence is difficult to follow. Consider restructuring it for clarity (e.g. “This makes the representation … complex and challenging.”).



L. 87, 204, 250, 273, Figure 1 caption: The manuscript frequently refers to a “macropore representation” or “macropore configuration”. Given that the approach does not explicitly represent macropore flow but rather mimics its effect on infiltration, I suggest using a more precise terminology (e.g. “conceptual macropore parameterization” or “macropore flow correction”) consistently throughout the manuscript.



L. 150: The statement that a unique residual unfrozen water content is assigned independently of soil temperature is somewhat unclear, as Eq. (1) explicitly includes temperature. I may be misinterpreting the formulation. However, it would be helpful to clarify whether this temperature dependence is averaged out to obtain a constant value for each soil type.

L. 214: The proposed parameterization represents a strong simplification of preferential flow processes. In particular it does not account for dynamic effects such as refreezing or the dependence of macropore flow on soil type and saturation. Some experimental studies even suggest that preferential flow can decrease at high saturation levels or vary significantly with soil structure. A short discussion (maybe later in the manuscript) of these limitations would strengthen the manuscript.



L. 311: The exclusion of stations based on NSE/KGE thresholds is understandable. However, it would be helpful to briefly explain the choice of the KGE threshold (-0.41) and its interpretation, as this value is not immediately intuitive for the readers.

L. 359: There appears to be a missing “Fr-” before the word “experiment”.

Figure3/Figure 5: The spatial distribution of excluded stations (red crosses) shows a pronounced clustering in specific regions rather than a random pattern. While the manuscript discusses potential causes for poor model performance in these areas, it remains unclear whether the station filtering systematically removes specific hydro-climatic or physiographic conditions from the evaluation. Please clarify how this affects the representativeness of the results.

L. 372: Please consider using more precise terminology such as “macropore parameterization” instead of implying an explicit representation.

L. 406: Please clarify whether this statement applies to all stations or only the majority of them.

L. 408: The underlying reasons for the observed behavior are not entirely clear. Additional explanation would improve the interpretation.

L. 414: Could the observed spatial variability in model performance be linked to differences in land use or soil properties? While some aspects are discussed, it might be helpful to further elaborate on the controlling factors (e.g. soil texture, land cover, or hydrological regime).

L. 418: The results indicate substantial performance degradations associated with the soil freezing configuration in some regions. While general limitations of the freezing scheme are discussed later, it would be helpful to more explicitly link these limitations to the observed regional patterns in model performance.

L. 421: The influence of station exclusion on the interpretation of PBIAS should be clarified, especially since different configurations are evaluated on different subsets of stations. This may affect the comparability of the reported bias between experiments. In addition, the use of the term “demonstrates” appears too strong in this context. Given the differences in station subsets and the indirect nature of the comparison, “suggests” may be more appropriate.

L. 439: While some explanations are provided, the underlying drivers of the observed spatial patterns are not fully clear. Could these patterns be more explicitly linked to dominant hydrological processes (e.g. snowmelt dynamics, flow partitioning, or model structural limitations)?

Figure 8: Adding a finer temporal reference (e.g. months or seasons) to the x-axis would improve readability and help distinguish between winter, snowmelt, and summer events.

Furthermore, differences between the model configurations appear to persist during summer periods, when soil freezing should not play a significant role. This may be related to the modified runoff generation scheme introduced together with the macropore parameterization. Please clarify whether such differences under unfrozen conditions are expected within the current model formulation.

Figure 9: Please clarify whether the improved agreement is related to specific freezing conditions during the selected period.

L. 610: The influence of soil type on the presented results could be discussed in more detail, as it may significantly affect the applicability and transferability of the parameterization. In particular, experimental studies suggest that the role of preferential flow under frozen conditions can vary with soil properties and saturation state. A brief discussion of this variability would strengthen the manuscript.

L. 615: The macropore configuration is implemented together with a modification of the runoff generation scheme by disabling the subgrid-scale interflow routine at the surface. This represents an additional structural change to the model. It would be helpful to more clearly separate and discuss the respective impacts of this modification and the macropore parameterization on the results.

General assessment

This manuscript addresses an important problem for cold-region hydrology and operational forecasting, namely the strong sensitivity of simulated streamflow to the representation of soil freezing. The study clearly shows that activating the soil-freezing configuration in SVS (Fr) strongly degrades hydrological performance, and that the Fr-MP configuration recovers a substantial part of that degradation. In that sense, the work has clear applied and operational value.

However, the manuscript currently overstates the physical meaning of the proposed approach and frames it as a macropore-process study. In my view, the paper does not demonstrate a macropore effect on streamflow, and the catchment-scale references cited in the Introduction do not establish macropores as the explanation for the weak runoff sensitivity to soil frost.

The manuscript would be much more convincing if it abandoned that framing and instead presented the proposed approach as what it is better described as: a conceptual correction of frozen-soil infiltration/runoff partitioning in an operational model, combining a relaxation of infiltration restriction under frozen conditions with a structural adjustment of runoff generation. More generally, I do not think the manuscript should adopt a macropore-based justification in the Introduction or present the parameterization in the Methods as if it were a process-based macropore representation. If the improvement brought by Fr-MP is to be related to real hydrological processes, this should instead be examined later in the Discussion, on the basis of the results and in a much more cautious and multifactorial way, including not only preferential flow, but also soil organic matter, peat and wetland hydrophysical properties, and broader hydropedological heterogeneity across the domain. The Title and Abstract should also be revised accordingly, so that they reflect this conceptual and operational framing rather than implying a process-based macropore study.

I therefore recommend major revision.

Major comments

1. The manuscript is framed too strongly as a macropore-process study

My main concern is that the manuscript is framed too ambitiously as a study of macropore effects on streamflow under frozen-soil conditions, whereas the results mainly show that the Fr-MP configuration pragmatically corrects a large part of the degradation introduced by Fr, without outperforming the noFr reference overall. The authors themselves show that Fr-MP improves Fr at many stations, but that noFr remains the best overall configuration. The Supplementary Material points in the same direction, since the additional hydrographs do not show a general superiority of Fr-MP over noFr.

This framing issue already appears in the Introduction. The manuscript builds a narrative according to which the discrepancy between strong infiltration restriction observed at small scale and the weak runoff sensitivity observed at plot or catchment scale can be explained, at least partly, by soil macropores. This interpretation is plausible as a broad conceptual hypothesis, but it is not directly established by the catchment-scale studies cited at that point. In particular, the two references used there to support the weak runoff sensitivity at catchment scale, Lindström et al. (2002) and Stähli (2017), do not directly demonstrate that this weak hydrological response is explained by macropores. Lindström et al. mainly conclude that no clear effect of soil frost on the timing and magnitude of runoff could be identified in their forested basin, and they emphasize that spring floods usually occurred when the soil was already unfrozen, together with an inverse relationship between frost depth and snow amount. Stähli (2017) likewise concludes that soil frost had no significant effect on winter runoff in the investigated pre-alpine catchments, but explicitly states that the study does not reveal why that effect is weak. He discusses several possible explanations, including shallow frost and strong spatial heterogeneity allowing concentrated infiltration, but without directly demonstrating that macropores are the dominant explanatory mechanism.

At the same time, the broader literature does support the idea that preferential or dual-domain flow may occur in frozen or structured soils at small scale, but it does so in a much more explicit and process-oriented way than the present manuscript. For example, Stähli et al. (1996) introduced a two-domain model approach for preferential flow in frozen soil, and Stähli et al. (1999) further showed, based on lysimeter experiments in two sandy soils, that infiltration into frozen soil depends on the full winter evolution of coupled heat and water dynamics and that key empirical parameters are not constant across soils and seasons. Demand et al. (2019) reported that connected biopores can promote bypass infiltration through shallow frozen layers under specific conditions, while Mohammed et al. (2018) reviewed this broader literature and emphasized the need for more explicit preferential-flow representations in frozen soils, including interacting macropore and matrix domains. More recently, Bauer et al. (2025), currently under review, supported the possibility of rapid bypass infiltration under frozen conditions in controlled experiments with artificial macropore networks, although referee comments on that preprint have also raised substantial concerns regarding over-interpretation of the results and the inferential nature of the proposed mechanisms.

Beyond the frozen-soil literature itself, other studies reinforce the same point. Smettem and Ross (1992) explicitly described a matrix-macropore dichotomy, showing that hydraulic conductivity can change by about an order of magnitude between slight tension and saturation, while localized preferred wetting may still occur even when matrix-based hydraulic properties correctly predict the onset of ponding. Baird (1997) likewise provided field evidence that macropores can be important for water and solute movement in fen peat and concluded that, if macropore flow is common, models based solely on Richards-type matrix flow may be inadequate and a matrix/macropore partitioning approach may be required. At the LSM scale, Vereecken et al. (2019) emphasized that infiltration remains difficult to upscale rigorously, that soil structural effects are still mostly neglected in land surface models, and that representing structure-related effects generally requires much more explicit treatment of soil heterogeneity and preferential pathways.

Taken together, the cited catchment studies do not establish a macropore-based explanation for the weak runoff response, whereas the process-oriented studies that do support preferential flow do so at much smaller scales and with much more explicit representations of matrix-macropore separation, soil structure, or preferential pathways. This reinforces my broader impression that the manuscript over-interprets its physical basis and overstates its “macropore effect” framing.

More fundamentally, the proposed formulation does not align with the process-based approaches commonly used to represent preferential flow: it neither represents organic-soil controls on hydraulic properties explicitly nor uses a dual-domain framework such as dual-porosity or dual-permeability. Instead, it applies a threshold-based removal of ice impedance together with a modification of runoff generation. It may therefore be useful as an operational correction, but it should be framed accordingly rather than as a process-based macropore study.

I therefore recommend that the manuscript be reframed explicitly, starting from the Introduction, as a conceptual correction of frozen-soil flow partitioning in an operational model, rather than as a process study of macropores.

2. Fr-MP combines multiple structural changes, so the gains cannot be attributed to “macropores”

Fr-MP does not merely alter a “macropore representation.” It combines at least two major structural changes: when W > W_MP, the ice-impedance effect is removed so that K_sat,v,fr = K_sat,v, and the surface-runoff rule is also changed so that runoff occurs only when the incoming vertical flux exceeds the local infiltration capacity, effectively suppressing the subgrid saturation-driven runoff mechanism at the surface. The discussion further states that the subgrid-scale interflow routine was disabled at the surface in Fr-MP.

As a result, the improvements cannot be attributed specifically to “macropores.” The paper shows that a combined modification of frozen-soil flow partitioning improves hydrographs relative to Fr. That is not the same as demonstrating a macropore effect. This point is made even stronger by the fact that infiltration itself is never shown explicitly as a diagnostic variable. The paper mostly shows changes in surface runoff, lateral flow, and soil drainage.

This concern also extends to the interpretation of the calibrated parameter α_MP. The calibration is performed over 1 September 2015 to 31 August 2018, with the first complete year used as spinup, whereas the main streamflow evaluation covers 2016–2021. Calibration and evaluation therefore overlap, which weakens claims about robustness and transferability of α_MP. In addition, α_MP is selected from a tested range extending from 0.5 to 0.99, and the optimal value retained by the authors (α_MP = 0.55) lies close to the most permissive end of that range. This suggests that the model mainly benefits from a strong relaxation of frozen-soil infiltration restriction, rather than from the identification of a physically robust macropore-activation threshold.

Finally, the statistical comparisons are not always based on identical station subsets. The three-experiment boxplot comparison is reported for 521 stations, whereas the manuscript also states that only 79% of stations are included in the Fr versus Fr-MP comparison. This does not invalidate the reported gains, but it does make the aggregate comparisons harder to interpret and weakens any strong claim that the calibrated α_MP has been shown to be spatially robust across the full domain.

3. The physical interpretation is limited by the soil formulation and the spatial heterogeneity of the domain

The physical interpretation of the results remains limited because the soil formulation used here is still essentially mineral-texture based, with hydraulic and freezing properties derived from variables such as X_sand, X_clay, W_sat, ψ_sat, and b. There is no explicit representation of organic soil horizons or organic-soil hydraulic properties, even though the study domain includes many wetlands and other cold-region environments where soil organic matter can strongly affect pore structure, storage, conductivity, and thermal behavior. Recent land-surface modeling work has shown that soil organic matter cannot be treated as a minor correction to mineral soil properties, and that physically consistent representation of organic matter is required to capture its effects on hydraulic and thermal processes. Decharme (2025) makes this point explicitly and shows that conventional LSM parameterizations based on mineral-soil assumptions or simplified SOC-based corrections can be physically inconsistent.

This limitation is particularly important in peat-rich or wetland environments. Liu and Lennartz (2019) show that hydraulic properties of peat soils depend strongly on pore structure, peat decomposition, and botanical composition, and that macroporosity and K_s follow different relationships with bulk density above and below a threshold near 0.2 g cm⁻³. Lennartz and Liu (2019) further emphasize that pristine peat is dominated by very high porosity and abundant macropores, whereas degradation reduces porosity and conductivity and may even turn highly degraded peat into a hydraulic barrier. Liu et al. (2020) then show that drainage progressively shifts peat pore structure over decadal to centennial timescales, reducing macroporosity, increasing small-pore volume, and decreasing K_s by roughly two orders of magnitude, with different trajectories under forest and agricultural land use. They also note that the different BD–K_s relationships cannot be explained by macroporosity alone, because pore connectivity and pore geometry also matter.

More broadly, recent work also shows that the spatial variability of K_s cannot be interpreted satisfactorily from soil texture and bulk density alone. Gupta et al. (2021) explicitly argue that texture-based pedotransfer approaches are limited because they neglect soil structure and pedogenic information, especially in vegetated and structured soils. Their results show that terrain, climate, and vegetation covariates affect spatial patterns of K_s, and that maps based on these environmental covariates better capture soil-formation-related spatial structure than maps based only on basic soil properties. At the same time, their spatial cross-validation performance remains modest, which underlines how difficult it is to interpret or predict K_s patterns robustly at large spatial scales.

In this context, I find the process-level interpretation in terms of “macropore effects” too broad. The model performs better in forested and more natural environments, and worse in several agricultural, urban, or hydrologically complex areas, but the current discussion does not convincingly disentangle what is due to frozen-soil physics, what is due to the mineral-soil formulation itself, and what may instead reflect wetland extent, peat degradation, drainage history, routing limitations, or other forms of hydropedological heterogeneity. This is especially important because the manuscript already acknowledges missing wetland representation and also includes non-explicit representations of tile drainage and ploughing, which further complicates any strict process interpretation in terms of macropores alone.

I therefore think the manuscript should be much more cautious in its spatial interpretation of the results. At present, the domain heterogeneity is too strong, and the soil formulation too simplified, to support a broad physical interpretation of the observed regional patterns in terms of macropores or preferential flow.

4. The failure of Fr is not diagnosed clearly enough

A central problem is that the paper still does not make it possible to diagnose clearly why the Fr configuration fails so strongly. To be fair, the manuscript does provide the main elements of the soil-freezing formulation itself, including the residual unfrozen water parameterization, the reduction of available pore space by ice, the impedance factor applied to hydraulic conductivity, and the threshold-based rule used in Fr-MP to suppress that impedance when W > W_MP ​. It also states that, in Fr-MP, surface runoff is only generated when the incoming vertical water flux exceeds the infiltration rate, thereby suppressing the subgrid saturation-excess runoff mechanism at the surface.

However, this remains insufficient for a clear mechanistic diagnosis of the Fr behavior, because the manuscript does not provide the full governing equations for the hydrological flux partitioning on which the interpretation actually relies. Surface runoff, lateral flow/interflow, and soil drainage are discussed extensively in the results and discussion, but their full formulations are not given explicitly in the manuscript itself. Instead, the reader is left with a conceptual diagram and a verbal description, while important structural details are only mentioned qualitatively, such as the disabling of the subgrid-scale interflow routine at the surface in Fr-MP and the strong near-surface anisotropy that promotes lateral flow.

As a result, it remains difficult to determine which component is primarily responsible for the severe degradation in Fr: the ice-impedance term itself, the runoff-generation scheme near the surface, the lateral-flow formulation, the anisotropy assumptions, or the interaction among these components. This matters because the paper interprets the improvements in terms of changes in surface runoff, lateral flow, and soil drainage, yet the reader cannot fully reconstruct how these fluxes are partitioned in each configuration. The fact that Fr-MP improves Fr while still generally underperforming noFr at the domain scale further reinforces the need for a more transparent diagnosis.

To make this diagnosis possible, I strongly recommend that the authors provide the missing governing equations, at least in the Supplementary Material. In particular, the manuscript should include the explicit formulations used for surface runoff generation, lateral runoff/interflow, and gravitational drainage, together with a compact summary of how these formulations are altered by Fr and by Fr-MP. That would make it much easier to understand which structural component of the frozen-soil configuration is actually responsible for the strong degradation and which component is effectively corrected in Fr-MP.

Minor comments

I do not have essential additional minor comments beyond those already raised by the other reviewers. However, I strongly suggest:

• not referring to this parameterization as a “macropore” representation, as it is more accurately described as a conceptual correction of frozen-soil infiltration/runoff partitioning in an operational model;
• adding the key equations of the model, especially those governing soil freezing, runoff generation, lateral flow, and drainage;
• clarifying whether the modified runoff-generation rule affects the model under unfrozen conditions;
• revising the abstract and conclusion so they more clearly reflect that Fr-MP mainly improves Fr, without outperforming noFr overall.

Additional references

Baird, A. J.: Field estimation of macropore functioning and surface hydraulic conductivity in a fen peat, Hydrol. Processes, 11, 287–295, https://doi.org/10.1002/(SICI)1099-1085(19970315)11:3<287::AID-HYP443>3.0.CO;2-L, 1997.

Decharme, B.: A process-based modeling of soil organic matter physical properties for land surface models – Part 1: Soil mixture theory, Geosci. Model Dev., 18, 9349–9384, https://doi.org/10.5194/gmd-18-9349-2025, 2025.

Gupta, S., Lehmann, P., Bonetti, S., Papritz, A., and Or, D.: Global prediction of soil saturated hydraulic conductivity using random forest in a Covariate-based GeoTransfer Function (CoGTF) framework, J. Adv. Model. Earth Syst., 13, e2020MS002242, https://doi.org/10.1029/2020MS002242, 2021.

Lennartz, B. and Liu, H.: Hydraulic functions of peat soils and ecosystem service, Front. Environ. Sci., 7, 92, https://doi.org/10.3389/fenvs.2019.00092, 2019.

Liu, H. and Lennartz, B.: Hydraulic properties of peat soils along a bulk density gradient—A meta study, Hydrol. Processes, 33, 101–114, https://doi.org/10.1002/hyp.13314, 2019.

Liu, H., Price, J., Rezanezhad, F., and Lennartz, B.: Centennial-scale shifts in hydrophysical properties of peat induced by drainage, Water Resour. Res., 56, e2020WR027538, https://doi.org/10.1029/2020WR027538, 2020.

Vereecken, H., Weihermüller, L., Assouline, S., Šimůnek, J., Verhoef, A., Herbst, M., Archer, N., Mohanty, B., Montzka, C., Vanderborght, J., Balsamo, G., Bechtold, M., Boone, A., Chadburn, S., Cuntz, M., Decharme, B., Ducharne, A., Ek, M., Garrigues, S., Goergen, K., Ingwersen, J., Kollet, S., Lawrence, D. M., Li, Q., Or, D., Swenson, S., de Vrese, P., Walko, R., Wu, Y., and Xue, Y.: Infiltration from the pedon to global grid scales: An overview and outlook for land surface modeling, Vadose Zone J., 18, 180191, https://doi.org/10.2136/vzj2018.10.0191, 2019.