General comments

This paper presents the development of a global 30-m seamless data cube by fusing Landsat and MODIS data, building on the authors' previous methods. This work is highly valuable for future applications requiring fine-resolution time series data. Despite the numerous algorithms developed for fusing Landsat and MODIS data or filling gaps in Landsat data, no global datasets generated by these technologies are currently available. However, the paper has several issues that need to be addressed.

Specific comments

1. Page 3, Line 70: For Landsat interpolation methods, there are techniques that do not require numerous clear-sky observations. For instance, the nearest similar pixel interpolator method only needs one clear-sky observation.
2. Introduction: Before the last paragraph, the authors should discuss the current research gap. Specifically, they should explain why a full-chain processing framework and such fused data are necessary. Additionally, they should outline the challenges users face when using current methods to produce data independently.
3. Page 11, Figure 3: When performing harmonization, was Landsat upscaled to the resolution of MODIS? Figure 3 suggests that both datasets have the same resolution. A linear transformation model was used; how do the authors address the issue when the linear model is not statistically significant?
4. Major Differences Between uROBOT and ROBOT Models: What are the primary differences between the uROBOT model and the previously developed ROBOT model by the authors?
5. Accuracy of the Time Series Model in Eq. 6: The accuracy of the time series model in Equation 6 could be affected by the time interval of the data. In cloudy regions, if the data is too sparse, is the result reliable? How do the authors address diverse changes when the data is sparse?
6. Eq. 7: The coefficients from MODIS are used for Landsat. This approach may be acceptable if the land is homogeneous, but it lacks a clear mechanism for complex landscapes. Landsat pixels have different temporal dependencies compared to coarse pixels. If this issue affects the reliability of the final product, the end users should be notified.
7. Eq. 9: The equation redistributes the residual to handle land cover changes, but the residual is at a coarse scale. How do the authors address changes at the fine pixel scale?
8. Page 15, Line 355: The three accuracy metrics mentioned cannot adequately assess the spatial context preserved by the fused images. Some spatial metrics, such as edge features (see examples in https://doi.org/10.1016/j.rse.2022.113002), should be presented.
9. Tables 4-6: Is the accuracy assessment conducted for each site? The tables only show the mean value. What is the range of these indices?
10. Figure 14: The RMSD values in Figure 14 are on a different scale compared to Tables 4-6.
11. Table 9: It would be beneficial to include a figure showing examples where only the Seamless Data Cube (SDC) can accurately classify the pixels, whereas other data cannot.

This submission aims at building a daily 30-m Landsat data record (2000-2022) based on Landsat 5,7,8,9 and MODIS for gap filling. Authors made sure that data are pre-processed (with AC and BRDF) and harmonized. The main pillars of the data is that they are daily at 30 m and span the 20+ years. While this dataset will certainly find a lot of applications and might be valuable to the community, and I applaud authors for  taking on this challenge, in my opinion, the description of this dataset is exaggerated (e.g., daily component) and authors do not provide evidences that this is actual daily data. You cannot use 8-(16-)-day data to prove that your dataset is actual daily because you should have available daily reference data.  I think this is one of the areas that authors did not work through and claim that your dataset resemble actual daily data is false. Maybe, it would have made sense to focus on regular 8-day composites. Maybe, the problem of generating daily data from 8-/16-day Landsat does not have a solution. Certainly, MODIS can help in certain regions, but I have huge doubts about its applicability globally given discrepancies in spatial resolution. Furthermore, even MODIS, majority of users use 8-day or 16-day composites given that MODIS acquires twice per day (Terra/Aqua). It's done because of clouds and because decrease in spatial resolution of the viewing angle (which increases with the cycle 1-15, and nadir only every 16 days). There are also some very dubious choices (e.g., moving the grid) . From CC comments, one can see obvious artifacts - it's understandable, from global product you can always find errors. But, as mentioned above, it's not real daily product, you did not provide evidences to claim it. I suggest authors to substantially re-work this and give a good thought on what real problem you are trying to solve and provide tangible solution that can be validated.



Issues:

1. "It is noteworthy that our adopted grid slightly deviates from the Sentinel-2 grid. Since the original Landsat coordinate system exhibits a half-pixel (15 meters) offset relative to the Sentinel-2 grid, we expanded and shifted the original MGRS grid by 15 meters in each direction to align with the Landsat coordinate system."

First, there is no shift; it's what is used for referencing a pixel value: center in Landsat and UL in S2 (in grid) (also what is used by default in GDAL). It was a poor decision to shift the grid. Whereas, if you selected MGRS as a coordinate grid, you should have re-projected Landsat (like HLS) into MGRS. Such shift can cause artifacts when comparing to the HLS data.



2. "Therefore, our approach aims at building multiple transformation models for each MGRS tile and each spectral band separately."

Building models per MGRS tile might introduce issues re spatial consistency. Did you check the impact of such an approach on the overlapping areas (between tiles)? How are consistent those temporally to reflect land cover changes?



3. I'm very skeptical about the applicability of MODIS gap-filling for Landsat on global scale. First, almost all existing approaches, including yours, does account for changes in spatial resolution with different angles. In reality, 500-m pixel can actually decrease up to 2 km one - see https://doi.org/10.1109/TGRS.2016.2604214

So, your assumption " the basic assumption of uROBOT is that the MODIS image 𝐶 can be accurately approximated by a linear combination of other similar MODIS images in the input timeseries data" is only valid under the condition that spatial resolution is invariant. That's not the case, especially in cloud-prone regions. Another issue that will not work in areas of small ag fields, e.g., Arica, SE Asia, etc.

I have seen multiple times examples when a bare ground field between two vegetative fields will be brightened in MODIS when resolution decreases. And used MODIS-based vegetation signal for Landsat will introduce huge errors. Therefore, a study must be conducted to explore how spatial resolution impacts restoration.



4. Section 5.3 must be re-written as it does not show limitations but rather what have not been done. There should be paragraphs re snow, re coastal regions and water, as AC algorithms do not work the best there (especially for snow as retrieval of aerosols is extremely difficult there). Furthermore, probably this product is not applicable to detecting rapid changes (daily) in land cover such as constructions as it will depend on actual acquisitions (and not blended daily). It will probably will not allow to detect "daily" burned fields or harvested fields because the signal will change in 2-3 days, or landslides, or iceberg movement, or fire propagation - anything that truly changes  every hour or day. Again, in reality your daily product will not allow (please, prove me wrong!) detection of these events (at daily basis) which require true daily data.