This is a well-structured paper that addresses a critical gap in Earth System Modeling (ESM) by proposing a novel stochastic approach to represent thermokarst lake dynamics at a regional scale. The model is conceptually sound, clearly explained, and the exploration of different dynamic regimes is a particular strength. The main limitations are openly acknowledged and are largely tied to the current state of observational data, which is a field-wide challenge, not a flaw specific to this study. Regarding the specific content of this manuscript, I have the following comments and suggestions.

1. The application of Geometric Brownian Motion (GBM) and Poisson processes to model lake area change and formation/drainage events is innovative and well-justified within the context of landscape-scale heterogeneity and stochasticity. This is a significant step beyond previous deterministic or analytical models.

2. The three idealized regimes (Complete Drainage, Oscillation, Stabilization) effectively demonstrate the model's behavioral range and provide a useful conceptual framework for understanding possible long-term trajectories of thermokarst landscapes. The links to real-world examples for each regime are well-made.

3. The current merging algorithm is identified correctly as a significant weakness. The assumption that merged lake area is the sum of the two original areas and that the new lake is perfectly circular is physically unrealistic and computationally expensive. It likely leads to a drastic overestimation of lake sizes and an underestimation of lake numbers.

4. The high percentage of unusable data points (30-33%) in the remote sensing dataset severely impacts the robustness of the parameter estimation. This uncertainty propagates through the observation-based simulations and limits the confidence in the derived parameters (λf, λd, μ, σ).

5. A crucial result is that the stochastic component (σ) dominates the deterministic drift (μ) in the observed period. This implies that, for the 2000-2020 period in this region, random environmental variability was a stronger driver of annual lake area changes than any clear climate-driven trend, explaining the lack of strong correlations.

6. The inability to find significant correlations between model parameters and climate variables (TDD, P) is a notable negative result. While honestly reported, it highlights the current impossibility of confidently projecting lake dynamics under climate change scenarios with this model, as intended in the abstract. This is a major constraint on its immediate application in ESMs.

7. The observation-based simulations project water area fractions increasing to over 50%, which is acknowledged as rare. This, combined with the high volatility, suggests the model parameters derived from 20 years of data may not be stable or representative of centennial-scale dynamics, potentially overestimating expansion.

8. The suggestion that the idealized simulations could be interpreted as spanning 10 ka with a 10-year time step is helpful for context, but the parameters would then be "per decade." This should be stated explicitly in the text to avoid confusion (e.g., in Table 1, add a note: "Parameters are per year; for a 10-year time step interpretation, values would be per decade").

9. The comparison with the van Huissteden et al. (2011) model is good. The explanation for the differing results (their reliance on a prescribed river network vs. your data-driven drainage rates) is plausible and highlights the advantage of your approach, but also its current data dependency.

10. The definition of "abrupt drainage" as a complete (>90% loss) and rapid event is clear. However, the discussion of results from other studies (Jones et al., 2011, 2020) that use different thresholds (e.g., >25% loss) is slightly confusing. A small table summarizing different study's definitions and converting their rates to a common framework would be helpful.

11. The axis titles of many figures are obscured, or the figures are not very clear, for example, in Figures 5, 6 and B1.

12. There are some spelling errors in the manuscript. For instance, the first reference should read "thermokarst lake" instead. The authors should carefully proofread the manuscript to avoid such errors.

Overall, this work is a valuable and necessary contribution to the field. It moves the conversation forward from deterministic single-lake models and analytical equilibrium solutions to a dynamic, stochastic, landscape-scale framework. While current data limitations restrict its predictive power, the model provides a powerful tool for hypothesis testing and a clear roadmap for what data is needed to eventually integrate these processes into climate projections.

The file is attached