RC1:

This manuscript addresses a long-standing challenge in the field of Land Data Assimilation: how to effectively assimilate Land Surface Temperature (LST) into offline Land Surface Models (LSMs), given that LST varies drastically over time and is predominantly driven by downward atmospheric radiation forcing. Based on the Common Land Model (CoLM) and the Local Ensemble Transform Kalman Filter (LETKF) algorithm, the authors propose a novel scheme to implicitly and jointly update surface soil temperature and soil moisture by utilizing the cross-correlations among variables.

Overall, this research topic holds strong scientific value and practical operational significance. The manuscript is logically structured with a rigorous experimental design, effectively overcoming the limitations of single-variable updates in traditional observation operators. It provides valuable insights for researchers in meteorology, hydrology, and Earth system modeling. However, to further enhance the scientific rigor of this study, there is still room for optimization and expansion. Specifically, the authors need to further justify the physical rationality of key parameter settings within the assimilation framework (such as the large-scale localization strategy, parameter perturbation magnitudes, and observation representativeness errors). Additionally, it is essential to demonstrate the comparative advantages of the LST observation operator, supplement the evaluation with independent in-situ validation, and provide an in-depth mechanistic interpretation of the degraded assimilation performance in specific regions. In light of these considerations, I recommend that the manuscript be accepted after Minor Revisions.

Reply:

We thank the reviewer for the insightful comments. Below is our point-by-point response to each concern.

RC1:

1. Localization Strategy in Land DA

The study adopts a distance-dependent decay localization strategy, with the localization radius D set to decrease linearly from 730 km at the equator to 146 km at the poles. While such a large-scale spatial smoothing strategy is common in atmospheric assimilation systems, the land surface exhibits extreme spatial heterogeneity (strongly influenced by topography, vegetation cover, soil texture, etc.). Forcing covariance propagation across different underlying surfaces (e.g., boundaries between lakes and deserts, or transition zones of different vegetation types) is highly prone to introducing spurious correlations. The authors are advised to thoroughly discuss the physical rationality of employing such a large spatial localization radius (up to 730 km) in land assimilation, and clarify whether matching constraints for underlying surface characteristics were considered during the covariance propagation process.

Reply:

The 730 km mentioned is the cutoff distance where the localization function drops to exactly zero [1]; the corresponding half-width is 365 km.

This localization scale is physically justified for two reasons:

1. It is considerably smaller than typical atmospheric system scales.

2. While larger than the minimum scale of land surface heterogeneity, it properly accounts for the spatial representation of aggregated satellite observations.

We have modified the text to add "cutoff" and explicitly clarified this rationale (L187-L195). Additionally, we noticed that the localization formula (Eq. (15)) in the original manuscript was mistakenly written in the form of an exponential function instead of a Gaspari–Cohn function, which has now been corrected in the revised manuscript.

RC1:

2. Model Parameter Perturbations

To maintain a reasonable ensemble spread, the authors introduced random perturbations to 14 tunable parameters in the CoLM model (see Table 1). Please provide further explanation regarding the physical justification for selecting these specific parameters and the rationality of setting their maximum perturbed magnitudes. Additionally, it is recommended to briefly discuss which of these 14 parameters contribute most significantly to the simulated variance of LST and the soil hydro-thermal cross-covariance.

Reply:

Table 1 in our manuscript summarizes the physical meanings for these 14 tunable parameters defined in the CoLM. For a more detailed explanation, please refer to the CoLM technical guide (http://globalchange.bnu.edu.cn/download/doc/CoLM/CoLM_Technical_Guide.pdf).

Regarding the perturbation magnitude, preliminary tests indicated that a 10% perturbation yielded an ensemble spread that was small relative to the observation error. We therefore increased the magnitude to 15% to generate a reasonable ensemble spread.

Owing to the highly nonlinear nature of land surface processes, any single one of these parameters can influence nearly all quantities in an LSM. Thus, it is not possible to isolate the sole contribution of each parameter. From a direct sensitivity perspective, LST is more sensitive to parameters related to surface energy balance and outgoing radiation (i.e., tuning factor, roughness lengths for soil and snow, drag coefficient, maximum transpiration, and critical temperature); and the hydro-thermal cross-covariance might be more directly governed by the parameters related to energy storage (i.e., maximum transpiration, critical temperature, wimp, irreducible water saturation of snow, restriction for minimum soil potential, pond depth, and fraction of model area with high water table).

RC1:

3. Lack of Independent In-situ Validation

The authors utilized ERA5-Land, GLDAS, and MERRA2 as evaluation benchmarks. Although these datasets offer the advantage of global coverage, they are inherently products driven by specific model parameterization schemes and atmospheric forcings, inevitably containing their own systematic biases. Please explain why independent in-situ observation data were not utilized for validation. If conditions permit, it is strongly recommended to supplement the study with site-level validation in representative regions to further enhance the objectivity and persuasiveness of the assimilation performance evaluation.

Reply:

Independent in-situ soil observations were not used for validation primarily due to the spatial scale mismatch and the lack of global soil observations. Therefore, spatially continuous datasets like ERA5-Land, GLDAS, and MERRA2 are more practical for large-scale benchmarking. Nevertheless, we agree that site-level validation would provide valuable insights, and we have explicitly discussed this limitation in the revised manuscript (L501-505).

RC1:

4. Observation Error Representation

Given the extremely high spatial resolution of MODIS LST, the authors used a simple area-weighted processing method and excluded cloud-contaminated pixels when aggregating it to the 0.5-degree model grids. It is recommended to further clarify how the observation error covariance matrix was explicitly specified in the LETKF update. When aggregating high-resolution observational data, did the observation error account not only for the retrieval error of MODIS itself but also reasonably incorporate the "representativeness error" caused by scale mismatch? This is crucial for the appropriate allocation of weights in the final analysis increments.

Reply:

As in other data assimilation studies, the observation errors are assumed to be independent, so the observation covariance is diagonal. The error variance (σ²) for each grid cell is a direct aggregation of the MODIS lst_uncertainty field, following the same quality control and aggregation procedures applied to the LST data.

We did not explicitly introduce a separate representativeness error term for scale mismatch. Preliminary tests indicated that the standard error of 0.01° LST with respect to the aggregated LST in a box of 0.5° × 0.5° is less than 5% of the aggregated retrieval uncertainty. This demonstrates that the representativeness error is negligible compared to the retrieval error.

RC1:

5. LST Operator

The accuracy of the observation operator itself plays a decisive role in the effectiveness of the assimilation. Since the method for calculating LST varies across different models, please explain the method used to calculate LST in CoLM and elaborate on its advantages or differences compared to observation operators in other models. This will help demonstrate the rationality of the adopted methodology.

Reply:

The simplest way to calculate radiometric LST is based on the net longwave radiation using the Stefan-Boltzmann relationship, LST_sim = (LW / σ)^(1/4). In CoLM, LW is constructed by combining thermal emissions from the sunlit canopy, shaded canopy, and the ground, while explicitly accounting for vegetation cover, canopy gap fraction, and surface emissivity. Eq. (2) is the formulation that the CoLM utilizes to calculate the LST, representing the realistic LST calculation  within CoLM's two-big-leaf framework [2]. We have added these details to the revised manuscript (L112-116).

RC1:

6. Degradation in Specific Regions

Figures 8 and 9 show that the simulation performance deteriorated in certain regions after LST assimilation (indicated by the red areas). Please provide an in-depth analysis of the physical mechanisms or algorithmic reasons leading to the degradation in these specific areas. Is this negative impact constrained by special local underlying surface conditions, limitations in the observation error settings, or a failure of the cross-correlations among physical variables? Please discuss whether this phenomenon is reasonable within the current assimilation framework.

Reply:

Degradation in soil water in western North America and northern Asia largely occurs in late spring or late autumn when freezing and thawing happen frequently. Assimilating the LST may alter the timing of thaw or freeze, which in turn causes the soil water to be discrepant from reference datasets.

A lake-land breeze-induced secondary circulation acts as a negative feedback mechanism that reduces the lake-land thermal contrast. However, in an offline LSM, the lack of this negative feedback can result in the degradation around the Great Lakes (largest in summer) because the assimilation updates only land grids while excluding water bodies.

The mechanism for northern Australia is already detailed in the manuscript. We have incorporated the analyses into the revised text (L424-428).

Reference:

[1] Gregory Gaspari and Stephen E. Cohn. “Construction of Correlation Functions in Two and Three Dimensions”. In: Quarterly Journal of the Royal Meteorological Society 125.554 (Jan.1999), pp. 723–757. issn: 0035-9009, 1477-870X. doi: 10.1002/qj.49712555417

[2] Yongjiu Dai, Robert E. Dickinson, and Ying-Ping Wang. “A Two-Big-Leaf Model for Canopy Temperature, Photosynthesis, and Stomatal Conductance”. In: Journal of Climate 17.12 (June 15, 2004), pp. 2281–2299. issn: 0894-8755, 1520-0442. doi: 10.1175/1520-0442(2004)017<2281: ATMFCT> 2.0.CO;2.

CC1:

This paper proposes a new way to assimilate satellite land surface temperature, specifically MODIS LST, into an offlineland surface model. The core idea is that LST should not update only temperature-related states, but should jointly update upper-layer soil temperature and soil moisture/water phase, using the cross-covariances inside an ensemble Kalman filtering framework. The authors implement this in CoLM with LETKF, apparently for the first time in this model, and run a global 0.5° experiment for 2001 with a 2-month free forecast afterward. The scientific motivation is reasonable: in offline LSMs, LST is hard to assimilate effectively because it changes rapidly and is strongly controlled by atmospheric forcing, so direct gains in LST itself tend to be short-lived. The manuscript argues that by storing the analysis signal in both energy and water states, the effect can persist longer and propagate vertically into deeper soil and snow states. Potentially worth publishing in HESS.

Reply:

We very much appreciate your valuable comments. In the following, we respond to each comment individually.

CC1:

The manuscript would benefit from a clearer statement of the exact novelty relative to Huang et al. (2008) and Chen et al. (2021). At present, the novelty is there, but it is somewhat buried in the introduction.

Reply:

The core distinction lies in the assimilation strategy. Our study proposes a joint assimilation framework that uses a single observation—Land Surface Temperature (LST)—to simultaneously update both surface soil temperature and soil moisture. By contrast, [1] assimilated LST to update only the soil temperature profile. Meanwhile, [2] assimilated LST for soil temperature and brightness temperature (BT) for soil moisture, but these act as independent updates rather than a joint assimilation.  We highlight the differences in the revised manuscript (L63-L66).

CC1:

Please clarify whether the 3-hour assimilation frequency is dictated by model cycling convenience or by the temporal handling of MODIS Terra observations, which are not naturally available every 3 hours globally. The observation-time treatment needs to be easier to follow.

Reply:

The 3-hour assimilation frequency represents a necessary balance between the physical characteristics of LST and model dynamics. Because LST is fast-varying, a narrow assimilation window is required to accurately capture its temporal changes. Conversely, the window cannot be too short, as the model requires sufficient spin-up time to dynamically adjust after each assimilation step. Therefore, a 3-hour cycle optimally satisfies both constraints.

For each 3-hour data assimilation cycle centered at 00:00, 03:00, ..., 21:00 UTC, only the satellite pixels whose actual scan times fall within this window are extracted and assimilated at that specific step. We have updated the manuscript to clarify this observation-time treatment (L210-L214).

CC1:

The threshold requiring at least 1500 raw MODIS observations per 0.5° grid cell seems high and may create regionally uneven sampling. Please justify this choice and show how much spatial coverage is lost because of it.

Reply:

The threshold of 1500 valid pixels (representing approximately 60% coverage) was chosen to ensure both spatial representativeness and high statistical confidence. A 0.5° grid cell in the CoLM model typically encompasses highly heterogeneous land cover and topography. Aggregating heavily clouded grids using an insufficient number of clear-sky pixels not only lacks statistical reliability but also inevitably introduces severe cloud-edge biases. Assimilating such unrepresentative observations would actively force the land surface model into an unrealistic physical state. Therefore, preserving physical consistency and statistical robustness must be prioritized over absolute spatial completeness.

To explicitly address the concern regarding coverage loss, we quantified the spatial gaps over a continuous one-month period (February 2001). We found that relaxing the threshold from 1500 down to 1000 valid pixels only marginally improves the spatial coverage of the observed grids by approximately 8.8%. This demonstrates that the data gaps are primarily driven by the physical reality of persistent cloud cover rather than the stringency of our algorithmic threshold. Consequently, compromising the representativeness and statistical standards would yield very limited gain in coverage while risking the introduction of significant assimilation errors.

We added the consideration for spatial representativeness in L221.

CC1:

The innovation rejection threshold of 8 K should be justified with either prior literature or a sensitivity test.

Reply:

We evaluated multiple threshold values (5 K to 9 K) to examine their impact on filtering the Observation-Minus-Background (OMB) innovations. It is worth noting that our assimilation window spans ± 1.5 hours. Given the pronounced diurnal cycle of land surface temperature (LST), this temporal span inherently introduces relatively large, yet physically valid, OMB deviations. Consequently, our tests indicated that a tighter threshold (e.g., 5 K) aggressively rejected a large portion of these valid observations, while a looser threshold (e.g., 9 K or higher) failed to adequately filter out gross errors and unrepresentative extremes. The 8 K threshold was empirically found to provide a good balance for the observation quality control in this study, effectively removing gross errors while safely retaining sufficient high-quality observations.

We have clarified the rationale behind this choice in the revised manuscript (L228-L229).

CC1:

The paper should report ensemble-spread diagnostics more explicitly, since ensemble collapse/spread deficiency is discussed as a central issue.

Reply:

To address this concern, we have added ensemble-spread diagnostics to the Discussion section (L488-L493). Taking the Northern Hemisphere as an example, the temperature spread in the first soil layer ranges from 0.3 K to 1.4 K. This relatively low spread is expected due to strong atmospheric forcing, which further justifies our approach of using LST to update both soil temperature and moisture. For soil moisture, the spread increases from approximately 0.5 kg/m² in the top layer to 10.0 kg/m² in the deepest layer, exhibiting depth-dependent and seasonal characteristics. Finally, as highlighted in the Conclusion, maintaining a reasonable ensemble spread in offline land surface models is fundamentally limited by atmospheric forcing. Overcoming this limitation requires a two-way coupled land-atmosphere assimilation framework.

CC1:

Please clarify whether perturbing the 14 tunable parameters is time-invariant per member or re-sampled through time. This matters for interpretation of model-error representation.

Reply:

The 14 tunable parameters are re-sampled at the beginning of each 3-hour assimilation cycle for each ensemble member rather than being time-invariant. We have added this clarification to ensure the model-error representation is accurately interpreted (L109-L110).

CC1:

I do strongly recommend that the authors cite Zafarmomen et al. (2024), “Assimilation of Sentinel-Based Leaf Area Index for Modeling Surface-Ground Water Interactions in Irrigation Districts”, as it is highly relevant to the present study.

Reply:

Reference [3] successfully assimilated high-resolution satellite-based Leaf Area Index into a coupled hydrological model to improve surface-groundwater simulations. This work is highly relevant to our study as it highlights the growing application of emerging remote sensing observations in hydrologic data assimilation. We have updated the Introduction (L29) accordingly.

Reference:

[1] C Huang, X Li, and L Lu. “Retrieving Soil Temperature Profile by Assimilating MODIS LST Products with Ensemble Kalman Filter”. In: Remote Sensing of Environment 112.4 (Apr. 15, 2008), pp. 1320–1336. issn: 00344257. doi: 10.1016/j.rse.2007.03.028

[2] Weijing Chen et al. “Retrieving Accurate Soil Moisture over the Tibetan Plateau Using Multisource Remote Sensing Data Assimilation with Simultaneous State and Parameter Estimations”. In: Journal of Hydrometeorology 22.10 (Oct. 1, 2021), pp. 2751–2766. issn: 1525-7541, 1525-755X. doi: 10.1175/JHM-D-20-0298.1

[3] Nima Zafarmomen et al. “Assimilation of Sentinel-Based Leaf Area Index for Modeling Surface-Ground Water Interactions in Irrigation Districts”. In: Water Resources Research 60.10 (Oct. 2024), e2023WR036080. issn: 0043-1397, 1944-7973. doi: 10.1029/2023WR036080

We sincerely thank the reviewer for the valuable time and insightful comments. Please find our detailed, point-by-point response addressing each of your concerns in the attached supplemental PDF.

We sincerely thank the reviewer for the valuable time and insightful comments. Please find our detailed, point-by-point response addressing each of your concerns in the attached supplemental PDF.