Widespread land surface cooling from paddy rice cultivation revealed by global satellite mapping

Domain: earth_science
Published: 2025-12-13
Editorial outcome: major_revision

Abstract
Paddy rice exacerbates climate warming through greenhouse gas emissions but also cools the land surface by enhancing evapotranspiration. While the former effect has received extensive attention, the biophysical cooling effect remains poorly quantified, partly due to the lack of high-quality global paddy rice data. Here, we address this gap by developing a universal rice mapping framework that integrates the strengths of phenology-based and curve-matching methods to construct the global, long-term rice dataset (GlobalRice500) with daily temporal and 500 m spatial resolution. Our analysis reveals that paddy fields annually reduce daytime land surface temperature by 0.21 (±0.0057)–0.27 (±0.0063) °C during the growing season compared to other croplands, with stronger cooling observed in larger fields and partial spillover to surrounding landscapes. These findings provide robust evidence of the surface cooling effect of paddy rice and call for a comprehensive evaluation of its role in climate regulation.

Introduction
Paddy rice is not only a cornerstone of global food security but also influences the global climate in at least two distinct ways. On one hand, rice cultivation accounts for nearly 30% of global agricultural methane emissions and ~48% of greenhouse gas emissions from croplands, contributing to global warming. On the other hand, paddy fields dissipate substantial heat through evapotranspiration, which helps lower land surface temperature (LST). However, previous studies have focused on the former (warming) effect but completely overlooked the latter (cooling) effect, potentially leading to a skewed understanding of the overall climate impacts of rice ecosystems. Elucidating the surface cooling effect of paddy fields is thus essential for a comprehensive assessment of their role in global climate regulation.

The global surface temperature effects of vegetation (e.g., forests and urban greening) and irrigation have been extensively studied. Although paddy fields consume approximately 30% of global irrigation water and represent a typical water–vegetation hybrid system, their cooling effect remains poorly quantified. The few existing studies are mostly limited to a relatively small geographic scope. One key limitation in addressing this knowledge gap is the lack of global rice distribution datasets with sufficient spatial and temporal resolution. As a comparison, forests typically exhibit long-term land cover continuity and relatively regular seasonal patterns. In contrast, paddy fields are subject to intensive human management, which often results in frequent land use changes and irregular seasonal dynamics. Such variability arises from practices such as crop rotation, temporary fallowing, or land conversion, as well as from the flexibility in planting schedules. In tropical regions, where water and heat are available year-round, rice can be cultivated at almost any time of the year. These characteristics make paddy fields especially challenging to monitor. Unlike forests, which can often be described using annual or multi-year averages, paddy fields require much higher temporal resolution to accurately capture the start, end, and duration of the rice growing season and to quantify their surface biophysical effect.

Considerable efforts have been made to map global rice area using satellite remote sensing since 2004 (Supplementary Table 1). However, to the best of our knowledge, none of the existing rice maps simultaneously satisfy the multiple requirements of global coverage, moderate spatial resolution, long temporal coverage and high temporal resolution. Most available datasets have a temporal resolution no higher than annual, with the exception of MIRCA2000 and GAEZ+_2015, which are monthly. Global-scale maps typically offer coarse spatial resolution (>5 km), and are often limited to a single year of temporal coverage, with the exception of GDHY. Higher spatial resolution maps generally have limited spatial coverage and are predominantly available for more developed regions, leaving data gaps in less-developed areas such as Africa, Central and West Asia.

Here, we construct a spatially and temporally comprehensive global rice dataset to assess the role of paddy rice in regulating land surface temperature. We address major bottlenecks in rice mapping caused by the variability in rice planting times and the deficiency of prior knowledge. We develop a universal rice mapping framework that identifies the start and end of the rice growing season on a daily basis, capable of recognizing paddy rice with arbitrary durations and intensities. This framework underpins the creation of GlobalRice500, a global-scale, long-term dataset of rice planting area with daily temporal and moderate spatial resolution. Finally, we analyze the cooling effect of paddy fields and characterize its spatial and temporal patterns, complemented by an uncertainty assessment to ensure the robustness of the results. Our findings offer a perspective and contribute to a more comprehensively understanding of the role of rice ecosystems in global climate regulation.

Results
The construction of GlobalRice500: long term mapping of global rice planting area
We first address the key bottlenecks in rice mapping by developing a universal rice mapping framework that integrates Moving Phenological Detection and Dynamic Time Warping (MPD_DTW, Fig. 1). MPD_DTW identifies the start and end of rice growing season on a daily basis, enabling the recognition of paddy rice with arbitrary growth durations and intensities. By leveraging the concept of “moving detection”, the MPD_DTW framework combines the dual strengths of both phenology-based and curve-matching methods. This approach has the advantage to go beyond rice mapping to offer inspiration for mapping other crops with distinct phenological phases, such as winter wheat characterized by a green-up stage.

We applied the MPD_DTW to analyze over 40 terabytes of satellite imagery spanning 2000–2023 and created GlobalRice500—the global-scale, long-term coverage (2001–2022), daily temporal resolution, 500m spatial resolution dataset of rice planting area (Fig. 2a). Since 2001, the global rice area has exhibited a fluctuating upward trend (Fig. 2b), with a net expansion of over 21 million hectares by 2022. Spatially, rice cultivation is primarily concentrated in Asia, with additional distribution found in Africa, North America, South America, Europe, and Oceania (Fig. 2c). From 2001 to 2022, Asia remained the world’s dominant rice-growing continent, accounting for >85% of the global rice area (Fig. 2c). Africa, now the second-largest rice-growing continent, recorded the fastest growth in rice area among all continents, more than doubling its area since the early twenty-first century (Fig. 2d). In contrast, the Americas experienced a notable decline (Fig. 2d). Driven by abundant heat and water availability, equatorial regions showed a higher propensity for rice cultivation (Fig. 2a). Aside from two leading rice-producing giants, India (ranked first) and China (ranked second), most of the other top countries with the largest rice areas were also located in tropical and subtropical regions (Fig. 3a).

For reference, we compiled existing local-scale high-resolution rice maps from diverse sources, including publicly available governmental datasets and rice maps released by academic researchers. The reference datasets cover all rice-growing continents and capture the diversity of global cultivation systems such as single-, double-, and multi-cropping (Supplementary Table 2). The overall accuracies typically exceed 90%, except in southeastern China, where it reached 88.85% (Fig. 3a and Supplementary Table 3). Moreover, comparision of country-level rice area between GlobalRice500 and FAO statistical data for 113 countries revealed a high degree of consistency, with an R-squared value of 0.98 (Fig. 3b). The spatial distribution of paddy rice derived by GlobalRice500 also aligned well with existing high spatial resolution rice maps (Fig. 3c–e).

The global cooling effect of paddy fields
The daily temporal resolution of GlobalRice500 enables precise delineation of rice growing seasons across diverse cropping systems worldwide, which is essential for quantifying the surface biophysical effect of paddy cultivation. This cooling effect by paddy rice is manifested in the reduction of land surface temperature (LST), which is measured as the difference in LST (ΔLST) between non-paddy croplands and paddy fields during the rice growing season (see “Methods”).

Our findings demonstrate that rice cultivation exerts a widespread cooling effect on the Earth’s surface. Globally, the annual mean daytime ΔLST between non-paddy croplands and paddy fields during 2003–2020 fluctuated between 0.21 (±0.0057) °C and 0.27 (±0.0063) °C (Fig. 4a). In contrast, the nighttime ΔLST was much smaller and nearly negligible—typically close to zero—indicating that the nighttime LST of paddy fields closely resembles that of surrounding non-paddy croplands (Figs. 4b and 5b). Uncertainty analyses of annual means for daytime and nighttime are reported in Supplementary Tables 4 and 5, respectively. During the rice growing season, daytime ΔLST exhibited an inverted V-shaped pattern with time, reaching its maximum value around Day of Season (DOS) 60–100 (Fig. 5c), largely corresponding to the peak growth stage.

We quantified the relationship between daytime ΔLST and paddy proportion (Methods) within 10-km grids and found that regions with higher paddy proportions exhibited stronger cooling effects (Fig. 5a). Regression analysis showed that ΔLST increased by approximately 0.0014 °C for every 1% increase in paddy proportion. Daytime ΔLST was further analyzed across different paddy field sizes, revealing that larger paddy fields tended to exhibit stronger cooling (Fig. 5b–d). Specifically, the maximum ΔLST during the rice growing season increased from 0.18 (±0.0032) °C at size 1 to 0.69 (±0.042) °C at size 11 (Fig. 5c).

To further investigate the spatial extent of the cooling influence, surrounding non-paddy areas were divided into concentric square-shaped rings according to their distances from the paddy fields, termed non-paddy ladders (NPLs) (Methods). Regardless of paddy field size, NPLs closer to paddy fields consistently showed lower ΔLST, which gradually increased and stabilized when moving further away from the paddy fields (Fig. 5d). This pattern indicates that the cooling effect of paddy rice extends beyond paddy field boundaries and spills over into surrounding non-paddy areas. Such spillover cooling suggests that paddy fields may contribute to local land surface temperature reduction, potentially mitigating urban heat island effects in adjacent cities. To mitigate the spillover effect of paddy field cooling, NPL was excluded in the calculation of ΔLST in Fig. 4 and Fig. 5a–c, with non-paddy areas represented by NPL (Methods).

Methods
Data sources and preprocessing
In this study, all available ~500m 8-day Moderate Resolution Imaging Spectroradiometer (MODIS) surface reflectance products from 2000 to 2023 were utilized to extract land surface information. Normalized Difference Vegetation Index (NDVI) and Land Surface Water Index (LSWI) were calculated from surface reflectance using the following formulas, and recorded their Days of Year (DOYs) in the meantime.

NDVI = (NIR - R) / (NIR + R) (1)

LSWI = (NIR - SWIR) / (NIR + SWIR) (2)

where R, NIR, SWIR represent the band of red (620–670nm), near infrared (841–876nm), short-wave infrared (2105–2135nm), respectively. NDVI could reflect crop biomass and chlorophyll content, since it’s correlated with leaf area index (LAI), green biomass, and percent green vegetation cover. LSWI is sensitive to the total amount of liquid water in vegetation and its soil background, so it has been widely used to identify flooding signals during rice transplanting period. Maximum synthesis of NDVIs from MOD09A1 and MYD09A1 was then performed each 8-day interval, in order to reduce the contamination of clouds and rain. Since the sensitivity of flooding signals is crucial for rice mapping, all the available LSWIs were retained. To acquire daily VIs, we applied a Savitzky-Golay filter to the NDVIs and a linear interpolation to the LSWIs.

To improve mapping efficiency and accuracy, Land Cover (LC) Products were added to mask out non-cropland area. We selected three LC products: IGBP (MCD12Q1 Land Cover Type 1) from 2001 to 2022, GlobeLand30 in 2010 and Crop Data Layer (CDL) in 2008–2022. GlobeLand30 and CDL were resampled into the format of 500m MODIS tiles, and pixels with cropland no less than 50% were reserved. Resampled CDL cropland were used as the cropland mask of the USA. For other area outside the USA, we chose the union of IGBP cropland and GlobeLand30 cultivated land as the cropland mask.

Rice maps around the world were collected for parameter selection and accuracy assessment (Supplementary Table 2), which may come from officially released data, previous studies, ground surveys, visual interpretation, etc. Rice maps were processed to the same format as 500m MODIS tiles, consistent with NDVIs and LSWIs. When processing a rice map with spatial resolution much higher than 500m (for example, 10m), we first scan the whole image with a L×L detecting window to select vast pixels, that means, fragmentary pixels which are not inside any detecting window would be excluded. L denotes the length of detecting window. Then, we resampled the scanned high-resolution image to 500m, labeled 500m pixels with rice proportion more than 50% as rice. For rice probability surface map in Australia, we used PhenoRice 50 algorithm to select highly reliable rice pixels. Moreover, for calibration and validation, we obtained country-level statistics from Food and Agriculture Organizations of the United Nations (FAO).

The global seamless 1km daily land surface temperature (LST) dataset in 2003–2020 was adopted to analyze the cooling effect of paddy fields, which contains daytime (13:30) LST and nighttime (01:30) LST. The Digital Elevation Model (DEM) dataset incorporated to improve the reliability of calculation of cooling effect. The DEM dataset we used was Global Version Elevation (GVE). Both the LST and DEM datasets were resampled to the 500m MODIS tile format for consistency with GlobalRice500.

Rice mapping framework
The MPD_DTW framework (Fig. 1) consists of two core modules: Moving Phenological Detection (MPD) and Dynamic Time Warping (DTW). The MPD module detects the start of the season (SOS) and end of season (EOS) of the rice growing season on a daily, pixel-by-pixel basis. SOS is identified based on the flooding criterion (Eq. 3).

LSWI + F > NDVI (3)

where F represents the flood factor. This process continues daily until LSWI + F exceeds NDVI, marking the onset of transplanting and thus defining the SOS. Following this, MPD identifies EOS when a local minimum of NDVI is detected, indicating a significant decline in biomass and greenness after harvest. EOS detection is further constrained such that

EOS - SOS ∈ [T1, T2] (4)

where T1 and T2 denote the lower and upper bounds of the locally defined rice growth duration. If no day satisfies the conditions of EOS within this expected time window, the module resumes the search for the next SOS.

Once an SOS and an EOS are established, the NDVI curve from this interval is analyzed through the DTW module. This module calculates the dissimilarity between the detected NDVI curve and a predefined standard rice NDVI curve derived from the collected reference rice maps. To accommodate the variability in the lengths of rice growing seasons, we employed the classic curve-alignment algorithm DTW to compute the dissimilarity. DTW nonlinearly aligns two time series of varying lengths to find the optimal warping path. The cost of the optimal alignment can be recursively computed by

D(Ai, Rj) = δ(ai, rj) + min {D(Ai-1, Rj), D(Ai, Rj-1), D(Ai-1, Rj-1)} (5)

where Ai = a1, ..., ai is the subsequence of the observed NDVI curve within the identified SOS–EOS period, and Rj = r1, ..., rj is the subsequence of a standard rice NDVI curve. The DTW dissimilarity is then given by the distance at the end of the alignment, D(A, R), which reflects the overall similarity in biomass and greenness dynamics between observed pixels and standard rice phenology.

If the DTW dissimilarity falls below a predefined threshold, current NDVI curve would be classified as a rice curve. And the detection of SOS will immediately skip to the identified EOS, ensuring no overlaps between different rice growing seasons. If the dissimilarity failed to satisfy the condition, the detection of SOS would continue at the day after current SOS.

Accuracy assessments of rice mapping
To evaluate mapping performance, two approaches were used. One is to precisely compute accuracies in reference regions where rice samples are present, the other is to compare mapped areas with statistical data across all subregions where official statistics are available.

Producer’s Accuracy (PA), User’s Accuracy (UA), Overall Accuracy (OA) and F1 score were used for accuracy assessment in reference regions. The corresponding formulas are as follows:

PA = TP / (TP + FN) (6)

UA = TP / (TP + FP) (7)

OA = (TP + TN) / (TP + FN + FP + TN) (8)

F1 = 2 × UA × PA / (UA + PA) (9)

where TP denotes true positives (i.e., rice correctly predicted as rice), TN stands for true negatives (i.e., non-rice correctly predicted as non-rice), FP represents false positives (i.e., non-rice incorrectly predicted as rice), and FN refers to the false negatives (i.e., rice incorrectly predicted as non-rice).

Country-level validation was performed against FAO statistics. For each country and year, we aggregated the pixel-level satellite-derived results into harvested area by summing the area of all pixels mapped as harvested and counting distinct harvest events within the calendar year (i.e., a pixel harvested k times contributes k × pixel area). This aggregation matches the spatial and temporal scales of the FAO data. Specifically, the relative errors (δ) between mapped areas and statistical areas were computed:

δ = (Area map,k,c - Area stat,k,c) / Area stat,k,c × 100% (10)

where k denotes year, c represents country, Area map means the harvest areas from mapping results and Area stat refers to the harvest areas from statistics.

Measuring cooling effect of paddy fields
We quantified the cooling effect of paddy fields by calculating the difference in LST between non-paddy cropland and paddy fields during the growing seasons (Eq. 11).

ΔLST = LST non-paddy - LST paddy (11)

where LST non-paddy and LST paddy represent the LST of non-paddy and paddy, respectively. Unless specified, LST refers to the mean value of LST during the rice growing season. In the calculation of the above mean, we excluded the 95th percentage to mitigate the impact of outliers.

To enhance comparability and reduce uncertainty, the non-paddy areas were selected under the following constraints: (1) non-paddy pixels were selected exclusively from cropland to mitigate disturbances from the underlying surfaces; (2) the paddy and non-paddy fields cannot be too distant from each other to ensure similar climatic conditions. We established a maximum allowable distance of 20 km between paddy and non-paddy fields; (3) elevation differences between paddy and non-paddy fields should be controlled to reduce temperature variations caused by changes in elevation. In our study, all elevation differences were less than 50m; (4) the number of non-paddy pixels need to be sufficiently large, reducing random errors. In our analysis, at least 10 non-paddy pixels were included in each comparison.

To analyze the relationship between ΔLST and paddy proportion, the paddy proportion refers to the percentage of paddy rice within a grid composed of 20 pixels × 20 pixels (approximately 10 km × 10 km). All 10 km × 10 km grids across the globe are binned based on paddy proportion, with each bin representing a 1% interval. For example, bin k represents all grids where paddy proportion in the range (k-1%, k]. The mean ΔLST during 2003–2020 of each bin is used as the representative ΔLST for that bin. Grids with a paddy proportion smaller than 0.2 were excluded, as their irradiance effect can be better observed when vegetation fraction is higher than 0.2.

To represent non-paddy areas at different locations, the non-paddy ladder (NPL) was defined. Specifically, NPL d refers to the set of non-paddy pixels within a square-shaped ring with a distance of d pixels from the edge of the paddy field (Supplementary Fig. 1a). NPL d1-d2 denotes the non-paddy pixels located within a square-shaped annular area between d1 and d2 pixels from the edge of the paddy field (Supplementary Fig. 1b). Please note that the distances described in the subscripts are measured in units of MODIS pixels (~500m).

In this study, if no specific NPL is indicated, the “non-paddy” in LST non-paddy default to NPL 10-30. NPL 10-30 was chosen instead of NPL 0-30 (excluding NRL 0-9). When examining the ΔLST across different non-paddy ladders (NPLs), non-paddy area in LST non-paddy is specified as a particular NPL d (Supplementary Fig. 1c). For example, In the x-axis of Figs. 2, 3 represents NPL 3, i.e., NPL 3-4.

Uncertainty analysis of ΔLST
For ΔLST estimates at 500m resolution, we quantified their uncertainty using a block bootstrap procedure, to account for both spatial and temporal autocorrelation. Each bootstrap sampling unit was defined as the combination of a spatial block (64×64 pixels, ~32×32 km) and a temporal unit (one calendar year). The ΔLST mean was recalculated for each bootstrap replicate. Based on 10,000 replicates, the standard errors of the mean (SEMs) and 95% confidence intervals (95% CIs) were derived from the bootstrap distributions. For the paddy-proportion-based means (Fig. 5a), the data had already been aggregated to 20×20 pixels; therefore, uncertainty was assessed directly using a Student’s t test. The “±” values reported throughout this paper represent the standard error of the mean (SEM) unless otherwise specified.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This study addresses an important yet understudied aspect of global climate dynamics related to paddy rice cultivation. The\nmanuscript is well-organized, and the methodological innovations introduced (GlobalRice500 using MPD_DTW) are\ncommendable. However, several critical points require further clarification and improvement:\nClassification Technique Concerns:\nThe authors utilize NDVI and LSWI indices within the MPD_DTW method for classifying rice paddies globally. Given the\nconsiderable variability in rice cultivation practices and environments worldwide, the robustness of relying solely on NDVI\nand LSWI needs deeper scrutiny. It remains uncertain whether the MPD_DTW framework can universally capture all the\nvariations accurately. The authors should explicitly discuss potential limitations or variations that may arise due to diverse\nagro-climatic conditions.\nValidation with FAO Data:\nIn Figure 1b, the authors present high accuracy compared to FAO statistics. It is critical to address the reliability and\nlimitations inherent to FAO datasets explicitly. FAO data typically represent aggregated regional or national-level data,\ncontrasting significantly with pixel-level satellite data. A clear logical explanation of how discrepancies between pixel-based\nsatellite observations and regional statistics are reconciled is necessary. Furthermore, areas exhibiting exceptionally high\nvalues in Fig. 1b suggest elevated uncertainty. Given that most rice paddies are concentrated in Asia, detailed validation for\nthese regions, specifically against more precise regional or local-scale datasets, would significantly enhance the study’s\nrobustness.\nStatistical Validation of GlobalRice500:\nA more thorough statistical validation and comparison of GlobalRice500 are warranted. Beyond FAO data, validation using\nmultiple independent datasets or statistical approaches should be undertaken to substantiate claims of accuracy and global\napplicability comprehensively.\nData Accessibility and Distribution:\nThe manuscript mentions that the GlobalRice500 dataset can be made available upon request. Given the dataset's potential\nscientific value, it is strongly recommended that the authors make it openly available through established repositories rather\nthan limiting distribution through the corresponding author’s requests.\nUncertainty and Statistical Analysis of ΔLST:\nThe analysis of ΔLST (delta Land Surface Temperature) is insightful, but the manuscript lacks a rigorous statistical\nvalidation and uncertainty quantification for ΔLST calculations. The authors should incorporate uncertainty analyses to\nstrengthen confidence in these findings.\nComparison of Cooling and Warming Effects:\nWhile the manuscript emphasizes the cooling effect of paddy fields (estimated at 0.21-0.27°C), it is essential to include a\ncomparative discussion with the warming effects due to greenhouse gas emissions from rice cultivation. An integrated\ndiscussion on the net climate effect (warming versus cooling) and its implications for global climate policy would significantly\nenrich the manuscript.", "reviewer_2": "This study uses a newly developed global dataset, GlobalRice500, to investigate the biophysical cooling effect of paddy rice\ncultivation on land surface temperature (LST). GlobalRice500 provides high-resolution (500 m), long-term (2001–2022) and\ndaily coverage of rice fields. The authors demonstrate that, compared to other croplands, paddy rice reduces daytime LST by\n0.21–0.27°C during the growing season, primarily due to differences in evapotranspiration. This cooling effect varies\nspatially and is stronger in regions with higher rice coverage, exhibiting spillover effects into adjacent non-paddy areas. The\nstudy emphasises the importance of incorporating biophysical cooling into climate models to counterbalance the recognised\nwarming effects of rice-derived greenhouse gas emissions.\nThe creation of GlobalRice500 fills a major gap in global rice mapping by combining high temporal resolution with global\nspatial coverage. However, I have reservations about the methodology underlying the statements made regarding the global\nland surface cooling from paddy rice cultivation. The following major concerns must be addressed before the paper is\nsuitable for publication.\nMajor concerns:\n1. The reported surface cooling effect attributed to paddy rice cultivation (ΔLST of 0.21–0.27 °C) is an order of magnitude\nsmaller than the known root mean square error (RMSE) of the MODIS-based LST product used in this study (approximately\n2.0–2.5 °C for daytime observations, according to Zhang et al., 2022). Although the manuscript aggregates data across\nmultiple years and spatial windows, it does not provide any formal uncertainty quantification to demonstrate that the cooling\nsignal is statistically robust and distinguishable from observational noise.\n2. Due to the small magnitude of the ΔLST signal in relation to the known sensor error, it is crucial to report confidence\nintervals or standard errors. A sensitivity analysis across regions with varying LST accuracy would further enhance the\ncredibility of the results. Without this, it remains unclear whether the observed cooling is a robust physical signal or falls\nwithin the bounds of measurement uncertainty.\n3. The attribution of cooling to rice cultivation relies on assumptions of climatic uniformity that may not account for important\nconfounding factors (e.g. albedo, soil moisture, irrigation infrastructure). As shown in Supplementary Tables 2 and 3,\nvalidation is heavily skewed towards Asia and North America. Given uneven validation coverage and reliance on visual\ninterpretation in some regions (e.g. Madagascar), a clearer justification of causality or more transparent discussion of\nlimitations is needed.\n4. The study emphasises that the cooling effect of rice has been overlooked in relation to its greenhouse gas emissions.\nHowever, the manuscript does not quantify or compare the radiative forcing effects of surface cooling versus GHG warming.\nWithout this information, there is a risk of misinterpreting rice cultivation as having net climate benefits, which is contrary to\ncurrent IPCC findings.\n5. The manuscript presents a 22-year rice map, but the validation appears to be limited to a few recent years in each region.\nThere is no evidence presented on the consistency of classification accuracy over time, which raises concerns about\npotential temporal drift, changes in sensor quality, cloud cover, and so on.\nMinor comments and suggestions:\n• The negligible nighttime ΔLST effect is interesting but underexplored. Provide maps of nighttime LST effects to support\nclaims of negligible impact of nocturnal transpiration and evaporation.\n• Consider including a simplified schematic of the MPD_DTW framework in the main text for clarity.\n• Define thresholds for phenological detection (e.g. \"flood factor F\") and DTW dissimilarity more transparently.\n• Sensitivity to grid size should be discussed. 500 m MODIS resolution may miss smallholder plots, particularly in Africa and\nSoutheast Asia.\n• No field-based validation is mentioned. Include more recent rice-mapping studies for benchmarking where possible (e.g.\n• Discuss data limitations: address potential biases in MODIS LST data; 1-km resolution LST data for fine-scale analysis of\n500-m rice maps.\n• Consider contextualising findings to establish its broader relevance: What is the magnitude of rice cooling compared to that\nof other land use effects (e.g. forests). What are the implications for agricultural policy such as water use trade-offs or crop\ndiversification?\n• Include uncertainty ranges for cooling estimates to reflect spatial variability and interannual fluctuations.\n• While GlobalRice500 covers 2001–2022, validation years are inconsistent (e.g. Madagascar: 2019; Australia: 2010). This\nraises questions about temporal representativeness.\n• Avoid overstatements without evidence such as \"natural air conditioners\"\n• Correct minor typos in the text such as \"sizzze\" in Fig. 2 caption.\n• Consider making the GlobalRice500 dataset publicly available."}, "author_response": "Response Letter to Reviewers’ Comments\nWe are grateful to the reviewers for their insightful comments and constructive\nsuggestions, which have significantly improved the quality and depth of our manuscript.\nIn this document, the point-by-point responses are shown in blue, and the content from\nthe manuscript and supplementary materials is shown in red. In the main manuscript and\nsupplementary materials, revisions are highlighted with red font.\nReviewer #1:\nThis study addresses an important yet understudied aspect of global climate dynamics\nrelated to paddy rice cultivation. The manuscript is well-organized, and the\nmethodological innovations introduced (GlobalRice500 using MPD_DTW) are\ncommendable. However, several critical points require further clarification and\nimprovement:\nClassification Technique Concerns:\nThe authors utilize NDVI and LSWI indices within the MPD_DTW method for\nclassifying rice paddies globally. Given the considerable variability in rice cultivation\npractices and environments worldwide, the robustness of relying solely on NDVI and\nLSWI needs deeper scrutiny. It remains uncertain whether the MPD_DTW framework\ncan universally capture all the variations accurately. The authors should explicitly discuss\npotential limitations or variations that may arise due to diverse agro-climatic conditions.\nResponse: We thank you for raising this important issue, which also provides valuable\ndirections for our future research.\nThe MPD_DTW framework is designed as a universal rice identification method\nadaptable across diverse agro-climatic regions and cropping systems. In the revised\nmanuscript, we have conducted multi-source and multi-scale validations of the\nGlobalRice500 dataset and added a new figure (Fig.3) to illustrate the results. These\ninclude pixel-level accuracy assessments against existing high-resolution rice maps, as\nwell as national-scale trend comparisons with statistical data. The validation dataset spans\nall rice-growing continents and encompasses various cultivation regimes, including\nsingle-, double-, and multiple-cropping systems (Fig.3). These results demonstrate that\nthe framework can robustly capture the large-scale spatial and temporal variations of\npaddy rice.\nMoreover, we have expanded the Discussion section to explicitly acknowledge and\ndiscuss the potential limitations of MPD_DTW framework.\nLine 530-536 (Figures):\nFig. 3 | Validations of GlobalRice500. a, Overall accuracies of GlobalRice500 across nine reference\nregions compared with satellite-derived rice maps. b, Comparison of GlobalRice500 and FAO\ncountry-level statistics for 113 countries (2001–2022), where each dot represents a country–year\nobservation. c-e, Visual comparison of rice spatial distribution between high-resolution rice maps and\nthe ~500 m GlobalRice500 for Myanmar (2019; 96.1°E, 16.6°N) (c), the USA (2022; 121.8°W, 39.1°N)\n(d) and China (2019; 130.6°E, 46.4°N) (e). High resolution rice maps include the 10 m NESEA-\nRice101, the 20 m SEA202 and the 30 m Cropland Data Layer (CDL)3.\nLine 157-164 (Discussion):\nWhile NDVI and LSWI effectively capture flooding and greening dynamics in many\npaddy rice systems, they may underperform in regions with atypical water-management\nregimes and in dry direct-seeded systems, where flooding signals are often weak or\nabsent. Moreover, although the 500 m spatial resolution of GlobalRice500 marks an\nadvance over most existing global rice maps (>5 km, Supplementary Table 1), mixed-\npixel effects cannot be entirely avoided, as a single pixel may encompass multiple land-\ncover types, leading to spectrally mixed signals and reduced classification accuracy.\nSmall, fragmented, and heterogeneous rice plots are particularly difficult to delineate\naccurately, whereas large, contiguous paddies can be mapped with greater confidence.\nValidation with FAO Data:\nIn Figure 1b, the authors present high accuracy compared to FAO statistics. It is critical\nto address the reliability and limitations inherent to FAO datasets explicitly. FAO data\ntypically represent aggregated regional or national-level data, contrasting significantly\nwith pixel-level satellite data. A clear logical explanation of how discrepancies between\npixel-based satellite observations and regional statistics are reconciled is necessary.\nFurthermore, areas exhibiting exceptionally high values in Fig. 1b suggest elevated\nuncertainty. Given that most rice paddies are concentrated in Asia, detailed validation for\nthese regions, specifically against more precise regional or local-scale datasets, would\nsignificantly enhance the study’s robustness.\nResponse: Thank you for this insightful comment.\nWe have now clarified how FAO statistics were used for validation. In our analysis, pixel-\nlevel satellite-derived maps were aggregated into annual country-level harvested areas by\nsumming the area of all pixels mapped as harvested within each country and year. To\nensure consistency with FAO’s reporting framework, we also accounted for multiple\nharvests within a calendar year—so a pixel harvested k times contributes k × pixel area.\nThis approach harmonizes the spatial and temporal scales of the two datasets prior to\ncomparison.\nWe fully acknowledge the differences between pixel-level satellite observations and\nnational-level statistics: FAO data provide aggregated regional or national totals, whereas\nsatellite observations capture fine-scale spatial heterogeneity. Accordingly, our\ncomparison with FAO statistics is not intended as pixel-level validation, but rather as a\nconsistency assessment of large-scale spatial and temporal patterns at the country level.\nIn the revised manuscript, we have\n(i) added explanations in the Methods section to clarify this procedure; and\nLine 312-315 (Methods):\nCountry-level validation was performed against FAO statistics. For each country and year,\nwe aggregated the pixel-level satellite-derived results into harvested area by summing the\narea of all pixels mapped as harvested and counting distinct harvest events within the\ncalendar year (i.e., a pixel harvested k times contributes k × pixel area). This aggregation\nmatches the spatial and temporal scales of the FAO data.\n(ii) expanded the Discussion section to address the reliability and limitations inherent of\nthe FAO dataset, highlighting how these complementary data sources together provide a\nmore comprehensive perspective on global rice distribution.\nLine 165-170 (Discussion):\nStatistical data (e.g., FAO) and satellite-derived observations (e.g., GlobalRice500)\npossess distinct advantages and limitations, yet their complementary strengths can be\nleveraged to provide a more comprehensive understanding of the global rice distribution.\nBuilding on this complementarity, our evaluation of GlobalRice500 integrates both\nperspectives: FAO statistics provide aggregated national or regional totals useful for\nbenchmarking large-scale, long-term consistency, whereas satellite-derived maps deliver\nfine-scale spatial details over specific periods, essential for characterizing local\nheterogeneity.\nFurthermore, we conducted more detailed validations of GlobalRice500 and\nincluded a new figure (Fig. 3) to illustrate the results. As the referee correctly pointed out,\n“most rice paddies are concentrated in Asia”; therefore, we placed particular emphasis on\nthis region when validating our dataset in the revised manuscript. Specifically, we\nselected four representative subregions characterized by distinct rice-cropping systems to\ncompute confusion matrices which are used for validation purpose: single-season\n(“Northeast China” and “Japan and Koreas”), double-season (“Southeast China”), and\nmulti-season (“Southeast Asia”). For detailed comparisons, we further validated\nGlobalRice500 against two high-resolution datasets (the 10 m NESEA-Rice101 and the\n20 m SEA202) in Myanmar, and with NESEA-Rice10 in Northeast China. Thanks to the\nreferee’s comment, we believe these additional validation efforts clearly enhanced the\ncredibility of our dataset.\nLine 530-536 (Figures):\nFig. 3 | Validations of GlobalRice500. a, Overall accuracies of GlobalRice500 across nine reference\nregions compared with satellite-derived rice maps. b, Comparison of GlobalRice500 and FAO\ncountry-level statistics for 113 countries (2001–2022), where each dot represents a country–year\nobservation. c-e, Visual comparison of rice spatial distribution between high-resolution rice maps and\nthe ~500 m GlobalRice500 for Myanmar (2019; 96.1°E, 16.6°N) (c), the USA (2022; 121.8°W, 39.1°N)\n(d) and China (2019; 130.6°E, 46.4°N) (e). High resolution rice maps include the 10 m NESEA-\nRice101, the 20 m SEA202 and the 30 m Cropland Data Layer (CDL)3.\nStatistical Validation of GlobalRice500:\nA more thorough statistical validation and comparison of GlobalRice500 are warranted.\nBeyond FAO data, validation using multiple independent datasets or statistical\napproaches should be undertaken to substantiate claims of accuracy and global\napplicability comprehensively.\nResponse: Thank you for highlighting this important point. In the revised manuscript, we\nhave added a new figure (Fig. 3) to present our validation results in detail. Overall, two\ncategories of datasets and three complementary validation approaches were employed to\ncomprehensively assess the accuracy and robustness of GlobalRice500.\nDatasets used for validation:\n(i) Statistical benchmarks to assess country-level annual consistency; and\n(ii) High-resolution rice maps to evaluate local spatial details.\nValidation approaches include:\n(i) Interannual agreement: Scatterplots comparing country-level annual estimates with\nstatistical benchmarks (Fig. 3b);\n(ii) Classification accuracy: Confusion matrices derived from comparisons with high-\nresolution maps, reporting overall accuracy, user’s and producer’s accuracies, and F1\nscores (Fig. 3a; Supplementary Table 3); and\n(iii) Local spatial detail: Side-by-side, zoomed-in comparisons in representative regions\nto illustrate spatial agreement (Fig. 3c-3e).\nSupplementary Table 3. Accuracies of reference regions\nRegion OA F1 UA PA\nnortheastern\n97.21% 90.51% 90.59% 90.43%\nChina\nsoutheastern\n88.85% 85.19% 85.63% 84.76%\nAsia China\nKoreas and\n93.45% 81.57% 81.85% 81.29%\nJapan\nSouthAsia 95.40% 89.82% 89.86% 89.78%\nOceania Australia 99.60% 79.27% 80.10% 78.46%\nAfrica Madagascar 95.32% 93.88% 93.88% 93.88%\nLombardia,\nEurope 94.42% 81.54% 81.75% 81.34%\nItaly\nNorth United States\n95.80% 80.73% 80.66% 80.80%\nAmerica of America\nSouth Rio Grande do\n94.34% 83.13% 83.26% 83.01%\nAmerica Sul, Brazil\n* OA, F1, UA and PA represent overall accuracy, F1 score, user accuracy and producer\naccuracy, respectively.\nData Accessibility and Distribution:\nThe manuscript mentions that the GlobalRice500 dataset can be made available upon\nrequest. Given the dataset's potential scientific value, it is strongly recommended that the\nauthors make it openly available through established repositories rather than limiting\ndistribution through the corresponding author’s requests.\nResponse: We appreciate your valuable suggestion. The GlobalRice500 dataset has been\nuploaded to Zenodo (DOI: 10.5281/zenodo.17460919). The dataset is currently under\nrestricted access, but it will be made publicly available immediately after the manuscript\nis accepted, to facilitate open access and broader scientific use.\nUncertainty and Statistical Analysis of ΔLST:\nThe analysis of ΔLST (delta Land Surface Temperature) is insightful, but the manuscript\nlacks a rigorous statistical validation and uncertainty quantification for ΔLST calculations.\nThe authors should incorporate uncertainty analyses to strengthen confidence in these\nfindings.\nResponse: Thank you for this valuable suggestion, and we apologize for the earlier\nomission. We have now incorporated a comprehensive uncertainty quantification for all\nΔLST estimates, with detailed procedures described in the Methods section.\nLine 355-364 (Methods):\nUncertainty analysis of ΔLST\nFor ΔLST estimates at 500 m resolution, we quantified their uncertainty using a\nblock bootstrap4,5 procedure, to account for both spatial and temporal autocorrelation.\nEach bootstrap sampling unit was defined as the combination of a spatial block (64×64\npixels, ~32×32 km) and a temporal unit (one calendar year). The ΔLST mean was\nrecalculated for each bootstrap replicate. Based on 10,000 replicates, the standard errors\nof the mean (SEMs) and 95% confidence intervals (95% CIs) were derived from the\nbootstrap distributions. For the paddy-proportion-based means (Fig. 5a), the data had\nalready been aggregated to 20 × 20 pixels; therefore, uncertainty was assessed directly\nusing a Student’s t-test. The “±” values reported throughout this paper represent the\nstandard error of the mean (SEM) unless otherwise specified.\n(1) Annual mean ΔLST (Fig. 4)\nUncertainty summaries for the annual daytime and nighttime means are reported in\nSupplementary Tables 4 and 5, including the standard errors of the mean (SEMs) and 95%\nconfidence intervals (CIs). And these have been referenced in the main text.\nLine 125-130 (Results):\nGlobally, the annual mean daytime ΔLST between non-paddy croplands and paddy fields\nduring 2003–2020 fluctuated between 0.21 (±0.0057) °C and 0.27 (±0.0063) °C (Fig.\n4a). In contrast, the nighttime ΔLST was much smaller and nearly negligible—typically\nclose to zero—indicating that the nighttime LST of paddy fields closely resembles that of\nsurrounding non-paddy croplands (Fig. 4b and Fig. 5b). Uncertainty analyses of annual\nmeans for daytime and nighttime are reported in Supplementary Tables 4 and 5,\nrespectively.\nSupplementary Table 4. Uncertainty analysis of global annual mean daytime ΔLST\nYear Mean (°C) SEM (°C) CI lower (°C) CI upper (°C)\n2003 0.2233 0.005774 0.2121 0.2348\n2004 0.2498 0.006438 0.2372 0.2624\n2005 0.2320 0.006214 0.2200 0.2444\n2006 0.2326 0.005756 0.2216 0.2438\n2007 0.2699 0.006332 0.2577 0.2822\n2008 0.2336 0.005727 0.2225 0.2450\n2009 0.2437 0.006050 0.2321 0.2557\n2010 0.2288 0.005855 0.2172 0.2402\n2011 0.2285 0.005929 0.2171 0.2404\n2012 0.2451 0.006217 0.2330 0.2572\n2013 0.2146 0.005712 0.2037 0.2260\n2014 0.2379 0.006058 0.2263 0.2498\n2015 0.2399 0.006035 0.2283 0.2517\n2016 0.2339 0.005870 0.2226 0.2456\n2017 0.2334 0.006113 0.2215 0.2455\n2018 0.2561 0.006076 0.2443 0.2681\n2019 0.2278 0.005988 0.2162 0.2396\n2020 0.2336 0.005908 0.2224 0.2455\n* SEM denotes Standard Error of the Mean. 95% CI lower and 95% CI upper represent\nthe lower and upper bounds of the 95% confidence interval.\nSupplementary Table 5. Uncertainty analysis of global annual mean nighttime\nΔLST\nYear Mean (°C) SEM (°C) CI lower (°C) CI upper (°C)\n2003 0.007687 0.001875 0.003956 0.01136\n2004 0.008449 0.002065 0.004474 0.01253\n2005 0.01808 0.001732 0.01470 0.02150\n2006 0.01311 0.001741 0.009729 0.01658\n2007 0.01298 0.002028 0.009019 0.01695\n2008 0.01819 0.001794 0.01457 0.02170\n2009 0.01251 0.001740 0.009115 0.01595\n2010 0.01664 0.001902 0.01292 0.02042\n2011 0.009387 0.001801 0.005846 0.01282\n2012 0.01318 0.001876 0.009485 0.01685\n2013 0.01103 0.001769 0.007612 0.01452\n2014 0.01390 0.001698 0.01057 0.01727\n2015 0.01575 0.001837 0.01208 0.01932\n2016 0.01306 0.001773 0.009581 0.01658\n2017 0.01167 0.001944 0.007949 0.01553\n2018 0.01073 0.001756 0.007292 0.01415\n2019 0.01405 0.001813 0.01050 0.01764\n2020 0.02231 0.001862 0.01863 0.02598\n* SEM denotes Standard Error of the Mean. CI lower and CI upper represent the lower\nand upper bounds of the 95% confidence interval.\n(2) Mean ΔLST across stratifications (Fig. 5):\nWe have revised Fig. 5 to visually display the uncertainty (SEMs or 95% CIs). The\nuncertainty quantification results for all panels are provided in the corresponding figure’s\nsource data file.\nLine 544-552 (Figures):\nFig. 5 | Multidimensional variations of ΔLST. a, Relationship between ΔLST and paddy proportion\n(Methods). Marker color indicates the standard error of the mean (SEM). b, Variations in ΔLST with\npaddy field size. Paddy size is represented by the side length (in pixels) of the paddy fields; for\nexample, a value of 5 denotes paddy fields measuring 5 × 5 pixels. Error bars indicate 95% confidence\nintervals (95% CIs). c, Seasonal variations of ΔLST across different Days of Season (DOS) for paddy\nfields of varying sizes, with curves smoothed using a five-day moving average. Numbers in\nparentheses indicate SEM. d, Spillover effect of paddy field cooling. The number k on the horizontal\naxis denotes the subscript of NPL (Methods). Different colors indicate paddy fields of varying sizes.\nk\nLines with markers show the means; shaded ribbons indicate the 95% CIs.\nComparison of Cooling and Warming Effects:\nWhile the manuscript emphasizes the cooling effect of paddy fields (estimated at 0.21-\n0.27°C), it is essential to include a comparative discussion with the warming effects due\nto greenhouse gas emissions from rice cultivation. An integrated discussion on the net\nclimate effect (warming versus cooling) and its implications for global climate policy\nwould significantly enrich the manuscript.\nResponse: Thank you for this valuable suggestion. In the revised Discussion, we have\nadded two new paragraphs comparing the cooling and warming effects of paddy fields\nand discussing their implications for global climate policy.\nNonetheless, it is important to note that these two effects are based on distinct physical\nconcepts and operate at different spatial and temporal scales. The warming effect, as\nreported in the IPCC Sixth Assessment Report (AR6), is expressed as changes in global\nmean surface air temperature (GSAT), resulting from the multi-decadal accumulation of\ngreenhouse gases in the atmosphere. In contrast, the cooling effect in our analysis is\nquantified using land surface temperature (LST), which represents biophysical response\nof surface energy fluxes that act locally.\nRigorous quantification and direct comparison between these effects are not yet feasible\ngiven their fundamentally different physical bases and spatio-temporal scales.\nNevertheless, we have provided indicative magnitudes to place these effects in context—\nglobal rice cultivation contributes roughly 0.03 °C to GSAT warming6,7, whereas our\nLST-based analysis indicates an annual biophysical cooling of 0.21–0.27°C (noting that\npaddy fields occupy only about 1% of the global land surface). We emphasize that future\nstudies should harmonize these metrics to better assess the integrated climate impacts of\nrice cultivation.\nLine 201-227 (Discussion):\nAs water-rich agroecosystems, paddy fields exert contrasting effects: abundant\nstanding water enhances evapotranspiration and induces local cooling, whereas persistent\nflooding creates anaerobic conditions that stimulate methanogenesis, thereby increasing\nmethane (CH₄) emissions8. According to the IPCC Sixth Assessment Report (AR6), the\nbest estimate of the total human-caused increase in global mean surface air temperature\n(GSAT) from 1850–1900 to 2010–2019 is 1.07 °C6. Global rice cultivation contributes\napproximately 2.5% to present-day anthropogenic warming7, corresponding to a GSAT\ncontribution of about 0.03 °C. These GSAT-based estimates reflect the global, multi-\ndecadal accumulation of atmospheric heat driven by greenhouse gases, whereas our LST-\nbased results quantify the global mean biophysical response of surface energy fluxes,\nprimarily driven by local evapotranspiration from paddy fields. Although our results\nreveal that paddy fields exhibit an annual mean cooling of 0.21–0.27 °C relative to other\ncroplands, it is important to note that they occupy only about 1% of the global land\nsurface9. Future studies should harmonize these temperature metrics to better quantify the\nintegrated climate impact of paddy cultivation.\nThese findings have broader implications for global climate policy. While\ngreenhouse gas mitigation remains the primary objective for reducing anthropogenic\nwarming, the biophysical cooling of paddy agroecosystems should also be recognized as\na local adaptation co-benefit10. Previous studies have shown that vegetation greening can\ninduce a net cooling effect on surface air temperature through enhanced\nevapotranspiration, partly offsetting greenhouse-gas-driven warming11. By quantifying\nthe biophysical cooling effects of paddy ecosystems, our study provides a complementary\nbasis for refining climate modelling and informing more balanced, climate-responsive\nagricultural strategies. At the field scale, such strategies could leverage the local cooling\nof paddy fields to buffer heat stress on crops and farm labor, thereby moderating hot-spot\ntemperatures during the growing season. In many regions, dryland crops in low-lying\nareas often suffer from waterlogging during heavy rainfall, whereas these same areas are\nnaturally suited for paddy cultivation, which can better utilize the standing water and\nmitigate heat stress. Where water availability allows, maintaining paddy in heat-prone\nsubregions could be considered to harness these microclimate benefits, though such\nmeasures must be carefully evaluated in light of greenhouse gas and water management\ntrade-offs.\nReviewer #2:\nThis study uses a newly developed global dataset, GlobalRice500, to investigate the\nbiophysical cooling effect of paddy rice cultivation on land surface temperature (LST).\nGlobalRice500 provides high-resolution (500 m), long-term (2001–2022) and daily\ncoverage of rice fields. The authors demonstrate that, compared to other croplands, paddy\nrice reduces daytime LST by 0.21–0.27°C during the growing season, primarily due to\ndifferences in evapotranspiration. This cooling effect varies spatially and is stronger in\nregions with higher rice coverage, exhibiting spillover effects into adjacent non-paddy\nareas. The study emphasises the importance of incorporating biophysical cooling into\nclimate models to counterbalance the recognised warming effects of rice-derived\ngreenhouse gas emissions.\nThe creation of GlobalRice500 fills a major gap in global rice mapping by combining\nhigh temporal resolution with global spatial coverage. However, I have reservations about\nthe methodology underlying the statements made regarding the global land surface\ncooling from paddy rice cultivation. The following major concerns must be addressed\nbefore the paper is suitable for publication.\nMajor concerns:\n1. The reported surface cooling effect attributed to paddy rice cultivation (ΔLST of 0.21–\n0.27 °C) is an order of magnitude smaller than the known root mean square error (RMSE)\nof the MODIS-based LST product used in this study (approximately 2.0–2.5 °C for\ndaytime observations, according to Zhang et al., 2022). Although the manuscript\naggregates data across multiple years and spatial windows, it does not provide any formal\nuncertainty quantification to demonstrate that the cooling signal is statistically robust and\ndistinguishable from observational noise.\nResponse: Thank you for raising this critical point. As the referee correctly pointed out,\nour analysis mitigates noise in the raw MODIS-based LST product through spatial and\ntemporal aggregation across multiple years and spatial windows. To formally address this,\nwe have now quantified the uncertainty of all aggregated ΔLST estimates using both\nstandard errors of the mean (SEMs) and 95% confidence intervals (CIs). We apologize\nfor the earlier omission and have revised the manuscript accordingly. Detailed procedures\nhave been added to the Methods section.\nLine 355-364 (Methods):\nUncertainty analysis of ΔLST\nFor ΔLST estimates at 500 m resolution, we quantified their uncertainty using a\nblock bootstrap4,5 procedure, to account for both spatial and temporal autocorrelation.\nEach bootstrap sampling unit was defined as the combination of a spatial block (64×64\npixels, ~32×32 km) and a temporal unit (one calendar year). The ΔLST mean was\nrecalculated for each bootstrap replicate. Based on 10,000 replicates, the standard errors\nof the mean (SEMs) and 95% confidence intervals (95% CIs) were derived from the\nbootstrap distributions. For the paddy-proportion-based means (Fig. 5a), the data had\nalready been aggregated to 20 × 20 pixels; therefore, uncertainty was assessed directly\nusing a Student’s t-test. The “±” values reported throughout this paper represent the\nstandard error of the mean (SEM) unless otherwise specified.\n(1) Annual mean ΔLST (Fig. 4)\nUncertainty summaries for the annual daytime and nighttime means are now\nreported in Supplementary Tables 4 and 5, including the SEMs and 95% CIs. We have\nalso added a corresponding description in the main text. The daytime ΔLST shows a\nconsistent and statistically significant cooling signal within narrow confidence intervals,\nconfirming that the observed daytime cooling signal is robust and distinguishable from\nobservational noise.\nLine 125-130 (Results):\nGlobally, the annual mean daytime ΔLST between non-paddy croplands and paddy fields\nduring 2003–2020 fluctuated between 0.21 (±0.0057) °C and 0.27 (±0.0063) °C (Fig.\n4a). In contrast, the nighttime ΔLST was much smaller and nearly negligible—typically\nclose to zero—indicating that the nighttime LST of paddy fields closely resembles that of\nsurrounding non-paddy croplands (Fig. 4b and Fig. 5b). Uncertainty analyses of annual\nmeans for daytime and nighttime are reported in Supplementary Tables 4 and 5,\nrespectively.\nSupplementary Table 4. Uncertainty analysis of global annual mean daytime ΔLST\nYear Mean (°C) SEM (°C) CI lower (°C) CI upper (°C)\n2003 0.2233 0.005774 0.2121 0.2348\n2004 0.2498 0.006438 0.2372 0.2624\n2005 0.2320 0.006214 0.2200 0.2444\n2006 0.2326 0.005756 0.2216 0.2438\n2007 0.2699 0.006332 0.2577 0.2822\n2008 0.2336 0.005727 0.2225 0.2450\n2009 0.2437 0.006050 0.2321 0.2557\n2010 0.2288 0.005855 0.2172 0.2402\n2011 0.2285 0.005929 0.2171 0.2404\n2012 0.2451 0.006217 0.2330 0.2572\n2013 0.2146 0.005712 0.2037 0.2260\n2014 0.2379 0.006058 0.2263 0.2498\n2015 0.2399 0.006035 0.2283 0.2517\n2016 0.2339 0.005870 0.2226 0.2456\n2017 0.2334 0.006113 0.2215 0.2455\n2018 0.2561 0.006076 0.2443 0.2681\n2019 0.2278 0.005988 0.2162 0.2396\n2020 0.2336 0.005908 0.2224 0.2455\n* SEM denotes Standard Error of the Mean. 95% CI lower and 95% CI upper represent\nthe lower and upper bounds of the 95% confidence interval.\nSupplementary Table 5. Uncertainty analysis of global annual mean nighttime\nΔLST\nYear Mean (°C) SEM (°C) CI lower (°C) CI upper (°C)\n2003 0.007687 0.001875 0.003956 0.01136\n2004 0.008449 0.002065 0.004474 0.01253\n2005 0.01808 0.001732 0.01470 0.02150\n2006 0.01311 0.001741 0.009729 0.01658\n2007 0.01298 0.002028 0.009019 0.01695\n2008 0.01819 0.001794 0.01457 0.02170\n2009 0.01251 0.001740 0.009115 0.01595\n2010 0.01664 0.001902 0.01292 0.02042\n2011 0.009387 0.001801 0.005846 0.01282\n2012 0.01318 0.001876 0.009485 0.01685\n2013 0.01103 0.001769 0.007612 0.01452\n2014 0.01390 0.001698 0.01057 0.01727\n2015 0.01575 0.001837 0.01208 0.01932\n2016 0.01306 0.001773 0.009581 0.01658\n2017 0.01167 0.001944 0.007949 0.01553\n2018 0.01073 0.001756 0.007292 0.01415\n2019 0.01405 0.001813 0.01050 0.01764\n2020 0.02231 0.001862 0.01863 0.02598\n* SEM denotes Standard Error of the Mean. CI lower and CI upper represent the lower\nand upper bounds of the 95% confidence interval.\n(2) mean ΔLST across stratifications (Fig. 5):\nWe have revised Fig. 5 to visually display the uncertainty (SEMs or 95% CIs). The\nuncertainty quantification results for all panels are provided in the corresponding figure’s\nsource data file.\nLine 544-552 (Figures):\nFig. 5 | Multidimensional variations of ΔLST. a, Relationship between ΔLST and paddy proportion\n(Methods). Marker color indicates the standard error of the mean (SEM). b, Variations in ΔLST with\npaddy field size. Paddy size is represented by the side length (in pixels) of the paddy fields; for\nexample, a value of 5 denotes paddy fields measuring 5 × 5 pixels. Error bars indicate 95% confidence\nintervals (95% CIs). c, Seasonal variations of ΔLST across different Days of Season (DOS) for paddy\nfields of varying sizes, with curves smoothed using a five-day moving average. Numbers in\nparentheses indicate SEM. d, Spillover effect of paddy field cooling. The number k on the horizontal\naxis denotes the subscript of NPL (Methods). Different colors indicate paddy fields of varying sizes.\nk\nLines with markers show the means; shaded ribbons indicate the 95% CIs.\n2. Due to the small magnitude of the ΔLST signal in relation to the known sensor error,\nit is crucial to report confidence intervals or standard errors. A sensitivity analysis across\nregions with varying LST accuracy would further enhance the credibility of the results.\nWithout this, it remains unclear whether the observed cooling is a robust physical signal\nor falls within the bounds of measurement uncertainty.\nResponse: Thank you for this constructive suggestion. In the revised manuscript, we have\n(i) quantified uncertainty by reporting 95% confidence intervals (CIs) and standard\nerrors (SEs) for all ΔLST estimates (Supplementary Tables 4–5 and Fig. 5), and\n(ii) added a sensitivity analysis across climate zones with varying LST accuracy\n(Supplementary Table 6) enhance the credibility of our result.\nAcross climate zones, the ΔLST estimates exhibit small SEs and narrow 95% CIs. All\ndaytime CIs lie above zero (cooling), which indicates that the observed cooling is a robust\nphysical signal.\nLine 187-191 (Discussion):\nGiven that the accuracy of satellite-derived LST retrievals varies by climate12, we further\nconducted a sensitivity analysis by stratifying ΔLST by Köppen–Geiger climate zones\n(2017)13. Zone-level means were summarized with 95% confidence intervals (CIs) and\nstandard errors (SEs) (Supplementary Table 6). The SEMs were consistently small and\nthe CIs were narrow across climate zones, indicating that, despite potential retrieval\nbiases, the observed cooling effect of paddy fields remains statistically robust.\nSupplementary Table 6. Uncertainty analysis across different climate zones\nRMSE of LST\nΔLST (°C)\nTime Climate zone (°C)\n25% 50% 75% Mean SEM CI lower CI upper\nEquatorial climate 1.97 1.96 2.03 0.1260 0.003566 0.1190 0.1330\nArid climate 1.40 1.42 1.64 0.4729 0.01196 0.4497 0.4962\nDaytime Warm temperate\n1.81 1.79 2.05 0.1372 0.003946 0.1296 0.1450\nclimate\nSnow climate 2.17 2.25 2.50 0.2592 0.01114 0.2380 0.2817\nEquatorial climate 1.24 1.21 1.33 0.01110 0.001400 0.008387 0.01386\nArid climate 1.03 1.19 1.30 0.01939 0.002917 0.01365 0.02520\nNighttime Warm temperate\n1.11 1.23 1.40 0.01773 0.001471 0.01485 0.02058\nclimate\nSnow climate 1.80 1.81 1.97 -0.02839 0.005069 -0.03829 -0.01855\n* The values of LST RMSE are quoted from the original publication12. 25%, 50%, and\n75% denote the exclusion rates.\n* SEM denotes Standard Error of the Mean. 95% CI lower and 95% CI upper represent\nthe lower and upper bounds of the 95% confidence interval.\n* Climate zones follow the Köppen–Geiger classification, 2017 version13.\n3. The attribution of cooling to rice cultivation relies on assumptions of climatic\nuniformity that may not account for important confounding factors (e.g. albedo, soil\nmoisture, irrigation infrastructure). As shown in Supplementary Tables 2 and 3, validation\nis heavily skewed towards Asia and North America. Given uneven validation coverage\nand reliance on visual interpretation in some regions (e.g. Madagascar), a clearer\njustification of causality or more transparent discussion of limitations is needed.\nResponse: We are grateful for your insightful comment. Our ΔLST calculation method\nfollows a commonly used standard strategy14,15,16. This approach compares the target\npixel(s) with adjacent pixel(s) under shared background conditions so that the observed\nLST contrast primarily reflects biophysical feedbacks rather than large-scale climatic\nvariability14. We acknowledge, however, that this design cannot fully eliminate the\ninfluence of factors such as albedo, soil moisture, or irrigation infrastructure, which may\nstill partially contribute to the observed LST differences. In the revised manuscript, we\nhave explicitly clarified these assumptions and limitations in the Discussion to ensure\ntransparent interpretation of our attribution.\nLine 178–181 (Discussion):\nTo ensure a comparable background, the non-paddy croplands serving as the control\ngroup were carefully selected under multiple constraints (Methods)14,15,16. Nevertheless,\nsome determinants of LST cannot be fully constrained in an observational design, and\nresidual differences in broadband albedo, soil moisture, and irrigation infrastructure may\nstill have influenced our estimates.\nRegarding the validation dataset, we note that the apparent spatial imbalance reflects both\ndata availability and the global distribution of rice agriculture. On one hand, Asia\naccounts for more than 85% of global rice harvested area and exhibits diverse cultivation\nsystems (e.g., single-, double-, and multi-cropping). Accordingly, we included multiple\nvalidation regions across Asia to capture this diversity and ensure that our assessment\nrepresents the full spectrum of global paddy management practices. On the other hand,\nthe United States provides comprehensive, publicly validated crop datasets, which are\nwidely used for benchmarking remote sensing products.\n4. The study emphasises that the cooling effect of rice has been overlooked in relation to\nits greenhouse gas emissions. However, the manuscript does not quantify or compare the\nradiative forcing effects of surface cooling versus GHG warming. Without this\ninformation, there is a risk of misinterpreting rice cultivation as having net climate\nbenefits, which is contrary to current IPCC findings.\nResponse: Thank you for this important comment. In the revised manuscript, we have\nexpanded the Discussion to compare the two effects and clarify their conceptual\ndifferences, thereby avoiding potential misinterpretation by readers.\nThe warming effect, as reported in the IPCC Sixth Assessment Report (AR6), is expressed\nas changes in global mean surface air temperature (GSAT), representing the multi-decadal\naccumulation of atmospheric heat driven by greenhouse gases. In contrast, the cooling\neffect in our analysis is quantified using land surface temperature (LST), which represents\nglobal mean biophysical response of surface energy fluxes that act locally through\nevapotranspiration.\nBecause these temperature metrics differ fundamentally in their physical meaning and\noperate at different spatial-temporal scales, a rigorous quantification or direct\ncomparability between these effects is not yet feasible. Nonetheless, we have provided\nindicative magnitudes to place these effects in context: global rice cultivation contributes\nabout 0.03 °C to GSAT warming, whereas our LST-based analysis indicates an annual\nbiophysical cooling of 0.21–0.27°C (noting that paddy fields occupy only ~1% of the\nglobal land surface).\nWe have explicitly cautioned readers that these magnitudes are not directly comparable\nand should not be interpreted as evidence of net cooling, but rather as complementary\ncomponents of paddy fields’ biogeochemical (radiative) and biophysical (non-radiative)\nclimate interactions.\nLine 201-213 (Discussion):\nAs water-rich agroecosystems, paddy fields exert contrasting effects: abundant\nstanding water enhances evapotranspiration and induces local cooling, whereas persistent\nflooding creates anaerobic conditions that stimulate methanogenesis, thereby increasing\nmethane (CH₄) emissions8. According to the IPCC Sixth Assessment Report (AR6), the\nbest estimate of the total human-caused increase in global mean surface air temperature\n(GSAT) from 1850–1900 to 2010–2019 is 1.07 °C6. Global rice cultivation contributes\napproximately 2.5% to present-day anthropogenic warming7, corresponding to a GSAT\ncontribution of about 0.03 °C. These GSAT-based estimates reflect the global, multi-\ndecadal accumulation of atmospheric heat driven by greenhouse gases, whereas our LST-\nbased results quantify the global mean biophysical response of surface energy fluxes,\nprimarily driven by local evapotranspiration from paddy fields. Although our results\nreveal that paddy fields exhibit an annual mean cooling of 0.21–0.27 °C relative to other\ncroplands, it is important to note that they occupy only about 1% of the global land\nsurface9. Future studies should harmonize these temperature metrics to better quantify the\nintegrated climate impact of paddy cultivation.\n5. The manuscript presents a 22-year rice map, but the validation appears to be limited to\na few recent years in each region. There is no evidence presented on the consistency of\nclassification accuracy over time, which raises concerns about potential temporal drift,\nchanges in sensor quality, cloud cover, and so on.\nResponse: We appreciate your valuable comment. To provide a clearer presentation of\nthe validation for GlobalRice500, we have added a new figure in the main text (Fig. 3).\nWe have made extensive efforts to collect and compile publicly available rice maps from\nmultiple sources, ensuring that all rice-growing continents are represented in the\nvalidation (Fig.3a). However, existing high-resolution datasets are generally limited in\nboth temporal and spatial coverage and are typically available only for a few years in each\nregion. As a result, our pixel-level validation is also constrained to those periods.\nThis limitation further underscores the motivation and contribution of our work—to\nproduce the first global, long-term dataset of rice area.\nWe also recognize the importance of evaluating the temporal consistency of classification\naccuracy. To address this, we aggregated our rice mapping results into annual harvested\nareas and compared them with the national-level statistics from the Food and Agriculture\nOrganization (FAO) for 2001–2022 to assess the consistency at the country scale (Fig.3b).\nThe strong agreement observed between our dataset and FAO statistics supports the\ntemporal robustness and reliability of our classification results.\nLine 530-536 (Figures):\nFig. 3 | Validations of GlobalRice500. a, Overall accuracies of GlobalRice500 across nine reference\nregions compared with satellite-derived rice maps. b, Comparison of GlobalRice500 and FAO\ncountry-level statistics for 113 countries (2001–2022), where each dot represents a country–year\nobservation. c-e, Visual comparison of rice spatial distribution between high-resolution rice maps and\nthe ~500 m GlobalRice500 for Myanmar (2019; 96.1°E, 16.6°N) (c), the USA (2022; 121.8°W, 39.1°N)\n(d) and China (2019; 130.6°E, 46.4°N) (e). High resolution rice maps include the 10 m NESEA-\nRice101, the 20 m SEA202 and the 30 m Cropland Data Layer (CDL)3.\nMinor comments and suggestions:\n• The negligible nighttime ΔLST effect is interesting but underexplored. Provide maps of\nnighttime LST effects to support claims of negligible impact of nocturnal transpiration\nand evaporation.\nResponse: Thank you for the suggestion. We have added the global map of nighttime\nΔLST (Fig. 4b). The map shows values close to zero across most regions, with only\nlocalized deviations, supporting our claims of negligible nighttime impact.\nLine 538-541 (Figures):\nFig. 4 | Global cooling effect of paddy fields. Geospatial distribution of mean daytime (a) and\nnighttime (b) ΔLST during 2003–2020, displayed as averages within 10 km × 10 km grids. ΔLST\nindicates the difference in LST between non-paddy croplands and paddy fields (non-paddy minus\npaddy) during the rice growing season.\n• Consider including a simplified schematic of the MPD_DTW framework in the main\ntext for clarity.\nResponse: We appreciate the recommendation and have incorporated a schematic of the\nMPD_DTW framework into the main text (Fig. 1) to clarify its underlying principles and\nworkflow.\nLine 511-521 (Figures):\nFig. 1| Schematic diagram of MPD_DTW. The diagram illustrates the workflow for identifying rice\ngrowing seasons based on daily time series data. The MPD module detects the start (SOS) and end\n(EOS) of the rice growing season (RGS) based on phenological knowledge. Thin vertical lines without\nsolid circles denote days that do not meet the identification conditions, whereas thick vertical lines\nwith solid circles indicate days that do. EOS is searched after the identified SOS (denoted as i in the\nfigure), within a search window from i + T1 to i + T2 (yellow block), which corresponds to the locally\ndefined rice growth duration. The NDVI curve within the detected SOS–EOS period (i to j) is then\nanalyzed by the DTW module to determine whether the time period corresponds to a RGS. The DTW\nmodule subsequently provides feedback to the MPD module, specifying the starting day for detecting\nthe next SOS. If the current curve is recognized as a RGS, the next SOS detection starts from the day\nafter the current EOS (j + 1), ensuring no temporal overlap between consecutive seasons. Otherwise,\nthe SOS detection restarts from the day following the current SOS (i + 1).\n• Define thresholds for phenological detection (e.g. \"flood factor F\") and DTW\ndissimilarity more transparently.\nResponse: Thanks for your suggestion. In the revised manuscript, we have defined the\nflood factor 𝐹 and the DTW dissimilarity more explicitly using equations.\nLine 272-275 (Methods):\nSOS is identified based on the flooding criterion (Eq. 3).\n𝐿𝑆𝑊𝐼 +𝐹 > 𝑁𝐷𝑉𝐼 , (3)\nwhere 𝐹 represents the flood factor.\nLine 285-293 (Methods):\nTo accommodate the variability in the lengths of rice growing seasons, we employed\nthe classic curve-alignment algorithm DTW to compute the dissimilarity. DTW\nnonlinearly aligns two time series of varying lengths to find the optimal warping path17,18.\nThe cost of the optimal alignment can be recursively computed by\n𝐷1𝐴 ,𝑅 4\n!#$ \"#$\n𝐷1𝐴 ,𝑅 4 = 𝛿1𝑎 ,𝑟4+𝑚𝑖𝑛< 𝐷1𝐴 ,𝑅 4 , (5)\n! \" ! \" ! \"#$\n𝐷1𝐴 ,𝑅 4\n!#$ \"\nwhere 𝐴 = 〈𝑎 ,…,𝑎 〉 is the subsequence of the NDVI curve of an observed pixel\n! $ !\nwithin the identified SOS–EOS period, and 𝑅 = 〈𝑟 ,…,𝑟〉 is the subsequence of a\n\" $ !\nstandard rice NDVI curve. The DTW dissimilarity is then given by the distance at the end\nof the alignment, 𝐷(𝐴,𝑅), which reflects the overall similarity in biomass and greenness\ndynamics between observed pixels and standard rice phenology.\n• Sensitivity to grid size should be discussed. 500 m MODIS resolution may miss\nsmallholder plots, particularly in Africa and Southeast Asia.\nResponse: We appreciate this comment and have clarified this limitation in the\nDiscussion: small, fragmented paddy plots are subject to higher detection uncertainty at\n500 m resolution.\nLine 160-164 (Discussion):\nMoreover, although the 500 m spatial resolution of GlobalRice500 marks an advance\nover most existing global rice maps (>5 km, Supplementary Table 1), mixed-pixel effects\ncannot be entirely avoided, as a single pixel may encompass multiple land-cover types,\nleading to spectrally mixed signals and reduced classification accuracy. Small,\nfragmented, and heterogeneous rice plots are particularly difficult to delineate accurately,\nwhereas large, contiguous paddies can be mapped with greater confidence.\n• No field-based validation is mentioned. Include more recent rice-mapping studies for\nbenchmarking where possible (e.g. 2023 Sentinel-1 rice maps at 20 m resolution).\nReference: Sun, C. et al. (2023). \"Twenty-meter annual paddy rice area map for mainland\nSoutheast Asia using Sentinel-1 synthetic-aperture-radar data.\" Earth System Science\nData, 15, 1501–1520.\nResponse: Thanks for the suggestion. We have strengthened the validation and added a\nnew figure (Fig. 3) in the main text. The figure includes benchmarking against recent\nhigh-resolution rice maps, including the 20 m Sentinel-1 product for mainland Southeast\nAsia (Sun et al., 2023) (Fig. 3c).\nLine 530-536 (Figures):\nFig. 3 | Validations of GlobalRice500. a, Overall accuracies of GlobalRice500 across nine reference\nregions compared with satellite-derived rice maps. b, Comparison of GlobalRice500 and FAO\ncountry-level statistics for 113 countries (2001–2022), where each dot represents a country–year\nobservation. c-e, Visual comparison of rice spatial distribution between high-resolution rice maps and\nthe ~500 m GlobalRice500 for Myanmar (2019; 96.1°E, 16.6°N) (c), the USA (2022; 121.8°W, 39.1°N)\n(d) and China (2019; 130.6°E, 46.4°N) (e). High resolution rice maps include the 10 m NESEA-\nRice101, the 20 m SEA202 and the 30 m Cropland Data Layer (CDL)3.\n• Discuss data limitations: address potential biases in MODIS LST data; 1-km resolution\nLST data for fine-scale analysis of 500-m rice maps.\nResponse: Thank you for this helpful suggestion. We have expanded the Discussion to\nacknowledge the data limitations, including potential biases in MODIS LST retrievals\nand scale mismatches between 1-km LST and 500-m rice maps\nLine 181-186 (Discussion):\nAdditionally, the LST dataset used predominantly represents clear-sky rather than\nall-weather conditions12. As it is derived from MODIS, retrievals rely on assumptions\nabout atmospheric profiles and surface emissivity, and residual biases may arise under\naerosol or thin-cloud contamination or over heterogeneous terrain19. Finally, combining\n1 km LST with 500 m rice maps inevitably entails a change of support; mixed-pixel and\nadjacency effects can introduce bias and attenuate field-scale LST gradients.\n• Consider contextualising findings to establish its broader relevance: What is the\nmagnitude of rice cooling compared to that of other land use effects (e.g. forests). What\nare the implications for agricultural policy such as water use trade-offs or crop\ndiversification?\nResponse: We appreciate your constructive suggestion. We have expanded the\nDiscussion to contextualize our findings within broader land-use and policy frameworks.\nSpecifically, we now\n(i) compare the magnitude of paddy-induced cooling with that of forest ecosystems; and\nLine 192-200 (Discussion):\nAmong different land-use and land-cover types, afforestation has been reported to\ninduce local cooling of 0.07–1.16 °C20. Satellite-based analyses further demonstrate that\nthe LST effects of forests vary by forest type and latitude—showing cooling of 0.27 ±\n0.03 °C in mid-latitudes, 0.50 ± 0.19 °C in southern high latitudes, and 2.41 ± 0.10 °C in\ntropical forests, but warming of 0.79 ± 0.03 °C in boreal forests15. By contrast, our results\nindicate that daytime mean cooling ranges from 0.13 ± 0.03 °C in smaller (~500 m × 500\nm) paddy fields to 0.43 ± 0.03 °C in larger (~5.5 km × 5.5 km) ones. It is worth noting\nthat the forest-related cooling reported above is measured relative to openlands, whereas\nour estimates are referenced to non-paddy croplands. Despite the smaller magnitude, the\npaddy-induced cooling represents a significant biophysical signal for an intensively\nmanaged and spatially fragmented agroecosystem.\n(ii) discuss the broader implications of paddy-induced cooling for agricultural planning\nand climate policy, including considerations of water-use trade-offs and adaptation co-\nbenefits.\nLine 214-227 (Discussion):\nThese findings have broader implications for global climate policy. While\ngreenhouse gas mitigation remains the primary objective for reducing anthropogenic\nwarming, the biophysical cooling of paddy agroecosystems should also be recognized as\na local adaptation co-benefit10. Previous studies have shown that vegetation greening can\ninduce a net cooling effect on surface air temperature through enhanced\nevapotranspiration, partly offsetting greenhouse-gas-driven warming11. By quantifying\nthe biophysical cooling effects of paddy ecosystems, our study provides a complementary\nbasis for refining climate modelling and informing more balanced, climate-responsive\nagricultural strategies. At the field scale, such strategies could leverage the local cooling\nof paddy fields to buffer heat stress on crops and farm labor, thereby moderating hot-spot\ntemperatures during the growing season. In many regions, dryland crops in low-lying\nareas often suffer from waterlogging during heavy rainfall, whereas these same areas are\nnaturally suited for paddy cultivation, which can better utilize the standing water and\nmitigate heat stress. Where water availability allows, maintaining paddy in heat-prone\nsubregions could be considered to harness these microclimate benefits, though such\nmeasures must be carefully evaluated in light of greenhouse gas and water management\ntrade-offs.\n• Include uncertainty ranges for cooling estimates to reflect spatial variability and\ninterannual fluctuations.\nResponse: Thank you for this helpful suggestion. We have now included uncertainty\nranges for all cooling estimates to capture both spatial variability and interannual\nfluctuations. Specifically, we report the standard errors and 95% confidence intervals for\n(i) annual ΔLST values (Supplementary Tables 4–5),\n(ii) stratified results across paddy proportion, paddy size, and non-paddy ladders (Fig. 5),\nand\n(iii) different climate zones (Supplementary Table 6).\n• While GlobalRice500 covers 2001–2022, validation years are inconsistent (e.g.\nMadagascar: 2019; Australia: 2010). This raises questions about temporal\nrepresentativeness.\nResponse: We appreciate your insightful comment. The variation in validation years\nacross regions primarily reflects the uneven temporal availability of high-resolution\nreference datasets. We have made extensive efforts to collect reference datasets, ensuring\nthat all major rice-growing continents are represented in the validation (Fig. 3a). However,\nexisting rice maps are typically produced for only a few years in each region.\nConsequently, our pixel-level validation is necessarily limited to those available periods.\nTo assess the temporal representativeness of GlobalRice500, we aggregated the mapped\nrice area into annual harvested areas and compared them with national-level statistics\nfrom the Food and Agriculture Organization (FAO) for 2001–2022 (Fig. 3b). The strong\ninterannual consistency between GlobalRice500 and FAO data supports the temporal\nrobustness and representativeness of GlobalRice500.\n• Avoid overstatements without evidence such as \"natural air conditioners\"\nResponse: Thank you for your suggestion. We have revised this sentence and carefully\nchecked the rest of the manuscript to avoid overstatements.\nLine 146-148 (Results):\nSuch spillover cooling suggests that paddy fields may contribute to local land surface\ntemperature reduction, potentially mitigating urban heat island effects in adjacent cities.\n• Correct minor typos in the text such as \"sizzze\" in Fig. 2 caption.\nResponse: Thank you for your reminding us, and we apologize for our oversight. We\nhave corrected “sizzze” to “size” in Fig. 2 caption, and have carefully proofread the entire\nmanuscript to eliminate other typographical errors.\n• Consider making the GlobalRice500 dataset publicly available.\nResponse: We appreciate your suggestion. The GlobalRice500 dataset has been uploaded\nto Zenodo (DOI: 10.5281/zenodo.17460919). The dataset is currently under\nrestricted access, but it will be made publicly available immediately after the manuscript\nis accepted, to facilitate open access and broader scientific use.\nReferences\n1. Han, J. et al. NESEA-Rice10: high-resolution annual paddy rice maps for Northeast\nand Southeast Asia from 2017 to 2019. Earth Syst. Sci. Data 13, 5969–5986 (2021).\n2. Sun, C. et al. Twenty-meter annual paddy rice area map for mainland Southeast Asia\nusing Sentinel-1 synthetic-aperture-radar data. Earth Syst. Sci. Data 15, 1501–1520\n(2023).\n3. USDA National Agricultural Statistics Service. Cropland Data Layer. USDA-NASS.\n4. Burke, M., Hsiang, S. M. & Miguel, E. Global non-linear effect of temperature on\neconomic production. Nature 527, 235–239 (2015).\n5. Moore, F. C., Baldos, U., Hertel, T. & Diaz, D. New science of climate change impacts\non agriculture implies higher social cost of carbon. Nat. Commun. 8, 1607 (2017).\n6. Intergovernmental Panel On Climate Change (Ipcc). Climate Change 2021 – The\nPhysical Science Basis: Working Group I Contribution to the Sixth Assessment Report\nof the Intergovernmental Panel on Climate Change. (Cambridge University Press,\n2023). doi:10.1017/9781009157896.\n7. Kritee, K. et al. High nitrous oxide fluxes from rice indicate the need to manage water\nfor both long- and short-term climate impacts. Proc. Natl. Acad. Sci. 115, 9720–9725\n(2018).\n8. Qian, H. et al. Greenhouse gas emissions and mitigation in rice agriculture. Nat. Rev.\nEarth Environ. 4, 716–732 (2023).\n9. FAO. FAOSTAT Statistical Database. Food and Agriculture Organization of the United\nNations.\n10. Bright, R. M. et al. Local temperature response to land cover and management\nchange driven by non-radiative processes. Nat. Clim. Change 7, 296–302 (2017).\n11. Zeng, Z. et al. Climate mitigation from vegetation biophysical feedbacks during the\npast three decades. Nat. Clim. Change 7, 432–436 (2017).\n12. Zhang, T., Zhou, Y., Zhu, Z., Li, X. & Asrar, G. R. A global seamless 1 km resolution\ndaily land surface temperature dataset (2003–2020). Earth Syst. Sci. Data 14, 651–664\n(2022).\n13. Rubel, F., Brugger, K., Haslinger, K. & Auer, I. The climate of the european alps:\nShift of very high resolution köppen-geiger climate zones 1800–2100. Meteorol. Z. 26,\n115–125 (2017).\n14. Li, Y. et al. Biophysical impacts of earth greening can substantially mitigate regional\nland surface temperature warming. Nat. Commun. 14, 121 (2023).\n15. Li, Y. et al. Local cooling and warming effects of forests based on satellite\nobservations. Nat. Commun. 6, 6603 (2015).\n16. Lee, X. et al. Observed increase in local cooling effect of deforestation at higher\nlatitudes. Nature 479, 384–387 (2011).\n17. Sakoe, H. & Chiba, S. Dynamic programming algorithm optimization for spoken\nword recognition. IEEE Trans. Acoust. Speech Signal Process. 26, 43–49 (1978).\n18. Petitjean, F., Inglada, J. & Gancarski, P. Satellite Image Time Series Analysis Under\nTime Warping. IEEE Trans. Geosci. Remote Sens. 50, 3081–3095 (2012).\n19. Wan, Z. New refinements and validation of the collection-6 MODIS land-surface\ntemperature/emissivity product. Remote Sens. Environ. 140, 36–45 (2014).\n20. Wang, H., Yue, C. & Luyssaert, S. Reconciling different approaches to quantifying\nland surface temperature impacts of afforestation using satellite observations.\nBiogeosciences 20, 75–92 (2023)."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The manuscript addresses the cooling effect of paddy rice cultivation well, and the use of GlobalRice500 with MPD_DTW is a strong contribution. The authors have clarified the limitations of NDVI and LSWI in different agro-climatic conditions and addressed concerns regarding data variability", "reviewer_2": "The revised manuscript provides a clear and reproducible description of the MPD_DTW framework, substantially improved global validation and well-structured comparisons with existing datasets that help establish the GlobalRice500 dataset. The key finding, i.e. consistent daytime surface cooling of 0.21–0.27 °C over paddy fields, is backed up by uncertainty analysis of the ΔLST estimates using block bootstrapping and reports confidence intervals that demonstrate the cooling effect despite the small magnitude relative to sensor error. The added sensitivity analysis across climate zones supports the consistency of results under varying LST retrieval accuracies. The discussion has been improved, now explicitly acknowledging residual confounding factors (e.g. albedo, irrigation and mixed-pixel effects) and contextualising the biophysical cooling within the broader greenhouse gas balance. While the revisions have substantially strengthened the manuscript, a few limitations remain that have only been partially addressed. Validation coverage remains uneven in some regions (e.g. Africa), and although confounding variables such as irrigation and albedo are discussed, they are not quantitatively controlled. The temporal stability of classification accuracy has not been demonstrated, nor have the implications of aggregation scale on ΔLST estimates been explored. While these limitations do not undermine the overall findings, they suggest that some claims of global representativeness and causal attribution should be interpreted with caution. These limitations are acknowledged."}, "author_response": "Response Letter to Reviewers’ Comments\nWe are grateful to the reviewers for their insightful comments and constructive\nsuggestions, which have significantly improved the quality and depth of our manuscript.\nIn this document, the point-by-point responses are shown in blue, and the content from\nthe manuscript and supplementary materials is shown in red. In the main manuscript,\nrevisions made in this round are highlighted with red font. Please note that only the\nchanges made in this round are highlighted; revisions from the earlier round are not\nmarked.\nReviewer #1 (Remarks to the Author):\nThe manuscript addresses the cooling effect of paddy rice cultivation well, and the use\nof GlobalRice500 with MPD_DTW is a strong contribution. The authors have clarified\nthe limitations of NDVI and LSWI in different agro-climatic conditions and addressed\nconcerns regarding data variability.\nResponse: We sincerely thank you for your encouraging comments on our study and\ngreatly appreciate your positive evaluation of our methodological contributions and\nclarifications.\nReviewer #2 (Remarks to the Author):\nThe revised manuscript provides a clear and reproducible description of the\nMPD_DTW framework, substantially improved global validation and well-structured\ncomparisons with existing datasets that help establish the GlobalRice500 dataset. The\nkey finding, i.e. consistent daytime surface cooling of 0.21–0.27 °C over paddy fields,\nis backed up by uncertainty analysis of the ΔLST estimates using block bootstrapping\nand reports confidence intervals that demonstrate the cooling effect despite the small\nmagnitude relative to sensor error. The added sensitivity analysis across climate zones\nsupports the consistency of results under varying LST retrieval accuracies. The\ndiscussion has been improved, now explicitly acknowledging residual confounding\nfactors (e.g. albedo, irrigation and mixed-pixel effects) and contextualising the\nbiophysical cooling within the broader greenhouse gas balance.\nWhile the revisions have substantially strengthened the manuscript, a few limitations\nremain that have only been partially addressed. Validation coverage remains uneven in\nsome regions (e.g. Africa), and although confounding variables such as irrigation and\nalbedo are discussed, they are not quantitatively controlled. The temporal stability of\nclassification accuracy has not been demonstrated, nor have the implications of\naggregation scale on ΔLST estimates been explored. While these limitations do not\nundermine the overall findings, they suggest that some claims of global\nrepresentativeness and causal attribution should be interpreted with caution. These\nlimitations are acknowledged.\nResponse: We sincerely appreciate the reviewer’s thorough and constructive comments,\nas well as the positive assessment of our methodological improvements, uncertainty\nanalysis, and expanded discussion.\nTo address the reviewer’s concerns, we have added further clarification in the\nDiscussion to help readers more objectively assess the spatial and temporal\nrepresentativeness of GlobalRice500.\nLine 176–181:\nHowever, our validation based on existing satellite-derived maps may exhibit spatial\nand temporal unevenness. Since Asia accounts for more than 85% of the global rice-\nharvested area and exhibits diverse cultivation systems, validation in this region is more\nextensive and detailed than in other regions. Moreover, existing high-resolution rice\nmaps are generally limited in temporal and/or spatial coverage (Supplementary Table\n1), resulting in an uneven and incomplete validation scope. This limitation further\nunderscores the motivation for, and contribution of, producing the first global, long-\nterm rice dataset, GlobalRice500.\nQuantitative control of confounding variables such as irrigation and albedo, as well as\nquantitative analysis of the scale dependence of ΔLST, are not yet feasible due to the\nlack of globally consistent, high-resolution datasets. We appreciate the reviewer’s\ninsightful comments on these issues, and future work will further investigate these\naspects as suitable data become available.\nReviewer #2 (Remarks on code availability):\nUnfortunately, the link provided did not work for me.\nError 404 was returned: This is not the web page you are looking for.\nResponse: The GitHub repository was previously set to private, which resulted in the\naccess error. It has now been made public and is fully accessible. We apologize for the\ninconvenience caused earlier."}]

Decision letter