Persistent impact of Fukushima decontamination on soil erosion and suspended sediment

In Fukushima, government-led decontamination reduced radiation risk and recovered 137Cs-contaminated soil, yet its long-term downstream impacts remain unclear. Here we provide the comprehensive decontamination impact assessment from 2013 to 2018 using governmental decontamination data, high-resolution satellite images and concurrent river monitoring results. We find that regional erosion potential intensified during decontamination (2013–2016) but decreased in the subsequent revegetation stage. Compared with 2013, suspended sediment at the 1-year-flood discharge increased by 237.1% in 2016. A mixing model suggests that the gradually increasing sediment from decontaminated regions caused a rapid particulate 137Cs decline, whereas no significant changes in downstream discharge-normalized 137Cs flux were observed after decontamination. Our findings demonstrate that upstream decontamination caused persistently excessive suspended sediment loads downstream, though with reduced 137Cs concentration, and that rapid vegetation recovery can shorten the duration of such unsustainable impacts. Future upstream remediation should thus consider pre-assessing local natural restoration and preparing appropriate revegetation measures in remediated regions for downstream sustainability. Mechanical soil decontamination is an important tool in remediating contaminated soils. Remediation efforts following the Fukushima Daiichi disaster increased soil erosion and downstream sediment loads that showed reduced 137Cs concentrations, but rapid revegetation quickly restored decontaminated landscapes.

R adioactive material leakage from nuclear industry activities or nuclear accidents poses a major threat to the environment and the economy 1,2 . Historically, widespread radioactive contamination has been observed several times, such as the Kyshtym accident (Soviet Union) in the 1940s-50s 3 , Windscale accident (England) in 1957 (ref. 4 ) and Chernobyl accident (Soviet Union) in 1986 (ref. 5 ). Long-term radiation exposure and radiophobia have increased the health risk and psychological pressure on the people of these regions, resulting in the abandonment of large areas rich in environmental resources and consequently in constrained sustainable human development 6,7 . As a key means for recovering contaminated regions, mechanical decontamination has been implemented in many legacy sites, including Hanford (United States) 8 and Chernobyl 9 . However, almost all attention has been directed to understanding in situ decontamination effects 10,11 and atmospheric particle resuspension issues 12 and little is known about if these perturbations would have secondary environmental impacts on their downstream catchments for the long-term.
On 11 March 2011, the most recent large-scale nuclear accident happened at the Fukushima Daiichi Nuclear Power Plant (FDNPP), Japan 13 . Over 2.7 PBq of fallout 137 Cs (half-life T 1/2 = 30.1 yr) from the FDNPP was deposited in the terrestrial environment, causing long-term radioactive contamination on large-scale neighbouring catchments 13,14 . To recover contaminated soil and decrease the exposure dose in Fukushima, the Japanese government evacuated the residents in 2011 and large-scale decontamination was implemented in contaminated villages 15,16 in 2012. Within a few years, dramatic land-cover changes occurred in the agricultural regions where 5 cm of the surface soil and vegetation were removed and replaced with uncontaminated soil (Fig. 1), with the subsequent natural restoration promoting revegetation in these decontaminated regions 11,13,16 .
The effectiveness of such intensive decontamination at reducing radiation exposure in situ is apparent, with air dose rates decreasing by 20-70% after decontamination 13 . However, as 137 Cs can firmly bind to clay minerals, it is transported along with suspended sediments (particulate 137 Cs) in the river system to the Pacific Ocean 17,18 . Comprehensively assessing the impacts of land-use changes in decontaminated regions on the downstream ecosystem is also necessary from the perspective of environmental sustainability. Moreover, recent studies have shown notable differences in the reduction in the particulate 137 Cs concentration over time across 30 rivers in Fukushima 19,20 , further underscoring the need to study land-use impacts on downstream particulate 137 Cs discharge into the ocean.
In addition to contaminant migration, land-use changes induced by strong perturbations (for example, decontamination) also alter land-ocean sediment transfer patterns 21,22 , which in turn affect elemental cycles 23 , biodiversity 24 and global climate 25 . Systematically assessing the perturbations' impacts on downstream river suspended sediments (SS) has thus become a joint goal of many related disciplines and a key part of developing science-based catchment management strategies [26][27][28][29] . However, owing to the limited availability of reliable and concurrent river monitoring data, how downstream river SS variations and land-use changes are linked remains unclear 21 .
With approximately 11.9% of the watershed area subjected to government-led decontamination between 2013 and early 2017, the Niida River Basin in Fukushima (Fig. 2a) 11,30 provides an excellent opportunity to examine the long-term impact on the dynamics of SS and particulate 137 Cs in river systems. Previous short-term river studies 19,20,[31][32][33][34] have suggested an increased river SS load and decreased particulate 137 Cs concentrations during decontamination. Several studies that analysed geochemical fingerprints in deposited sediments have implied a significant contribution of sediment sources from the decontaminated regions to the river [35][36][37] . Yet, given that threats to sustainable catchment management from anthropogenic perturbations are often long-term 28  Haramachi catchment, respectively. c, Ordered decontaminated area for agricultural land (including paddy field and cropland), residential land and grassland, by year. The digital elevation model (DEM) data were from the Geospatial Information Authority of Japan. The 137 Cs inventory data were obtained from Kato et al. 14 . The information for the scheduled decontamination regions (paper maps) was obtained from the Ministry of the Environment of Japan. and reliable data on quantified land cover and continuous river records are required to explore the effect of decontamination on river SS and particulate 137 Cs discharge.
Here we provide a comprehensive assessment of the impacts of land-use changes in decontaminated regions on river SS and particulate 137 Cs dynamics, as well as the downstream discharge. We mapped the evolution of decontaminated region boundaries using governmental decontamination documents (Fig. 2b), photographed the land cover in the decontaminated regions using drones and quantified the land-use changes using the normalized difference vegetation index (NDVI) at 10 m spatial resolution. Meanwhile, we conducted a long-term field investigation (Methods and Supplementary Table 1) spanning the decontamination (2013-2016) and natural restoration (2017-2018) stages to continuously record the fluctuations in water discharge and turbidity (10 min temporal resolution) and particulate 137 Cs concentrations, both upstream and downstream. Combining the above quantitative data, we systematically reveal that long-term land-use changes in upstream decontaminated regions greatly affect sediment and 137 Cs discharge from downstream river systems into the Pacific Ocean.

land-cover changes in decontaminated regions
Regional decontamination was accomplished in March 2017, spanning over 22.9% and 11.9% of the upstream (Notegami) and downstream (Haramachi) watershed areas, respectively (Fig. 2c). In 2014, agricultural land (18.02 km 2 ) was one of the major land uses in the regions where decontamination was ordered, with a significant increase of over 720% compared with that in 2012. Conversely, the changes in the ordered grassland (1.93 km 2 ) were minimal, with an approximately 26% increase between 2012 and 2014. Given that overall land-cover changes were more pronounced in agricultural lands than in grasslands or residential lands, large-scale agricultural land decontamination may severely alter landscape erodibility and consequently the sediment supply. Moreover, the proximity of the decontaminated agricultural land to rivers increases sediment transport from terrestrial environments.
Drone photographs ( Fig. 3a and Supplementary Fig. 1) showed significant land-cover changes in the upstream decontaminated regions. For instance, the agricultural land at the Hiso site (D3, Fig. 3a) was almost bare in August 2016 when the decontamination was completed. However, natural restoration caused the recovery of vegetation during the post-decontamination stage (August 2018). Considering the decontamination sequence and seasonal dependence of the plant growth cycle 38 , drastic spatiotemporal land-cover changes in the decontamination regions are conceivable.
To quantitatively estimate and compare land-cover changes in the decontaminated regions, we generated NDVI maps in these regions based on the available satellite images from Sentinel 2 and the Moderate Resolution Imaging Spectroradiometer (MODIS) between 2011 and 2018. The comparison between realistic scenarios and Sentinel 2-based NDVI maps (spatial resolution: 10 m) from drone observation sites ( Fig. 3a and Supplementary  Fig. 1) confirmed the feasibility of NDVI images in distinguishing bare land and agricultural land with vegetation cover. To improve the temporal-spatial resolution of the NDVI dataset, we used the enhanced spatial and temporal adaptive reflectance fusion model (ESTARFM) 39 to fuse Sentinel 2 and MODIS maps. Subsequently, we used interpolation to link these newly generated NDVI data (spatial resolution: 10 m) and plotted daily NDVI variations for all decontaminated regions between 2011 and 2018 (Fig. 3b).
In the daily NDVI variation curve, the peak NDVI during 2013-2014 was similar to the pre-decontamination stage (2012) but decreased by approximately 10% in 2016. After decontamination, the peak value presented an increasing trend under the influence of vegetation recovery. However, as the government lifted the evacuation zone after 2017 and allowed the residents to return, vegetation was again removed from some areas planned for agricultural activities in 2018 ( Supplementary Fig. 1b). Further analyses of the NDVI variations in the decontaminated regions scheduled in 2012 (Fig. 3c), 2013 (Fig. 3d) and 2014 (Fig. 3e) showed that the NDVI peaks were decreased by approximately 12%, 11% and 15%, respectively, within 2-3 years after decontamination was ordered, thereby providing unambiguous evidence for decreasing vegetation land cover caused by decontamination.
We converted all NDVI maps derived from ESTARFM-images to C × P (cover management and support practice factors, respectively) maps using empirical models 40,41 . We then estimated erosion potential (that is, K × LS × C × P in the revised Universal Soil Loss Equation (RUSLE); Methods) maps using the LS (slope length and slope steepness factors, respectively) map ( Supplementary Fig. 2) and K (soil erodibility) factor in the decontaminated regions. We found that the slopes of decontaminated regions were generally similar for each decontamination-ordered year ( Supplementary  Fig. 2), suggesting that the erosion potential was consistent with the NDVI in the decontaminated regions. Therefore, we also estimated the daily variation curve of the erosion potential using the mean CP, LS and K factors in the decontaminated regions.
Here we show the ESTARFM-based NDVI (Fig. 3f) and erosion potential (Fig. 3g) during the summer season for each year (specific periods in Supplementary Table 2). NDVI showed a decreasing trend from 2013 to 2016, while the erosion potential peaked in 2016, representing approximately 98% and 52% increases over the pre-decontamination (2011) and natural restoration (2018) stages, respectively. Combining the corresponding NDVI ( Supplementary Fig. 3) and erosion potential maps ( Supplementary  Fig. 4), significant changes in the spatial differences between the land cover and erosion potential during decontamination were also observed.

response of river SS to land-cover changes
The downstream river SS load (L; Fig. 4a) exhibited a strong correlation with water discharge (Q) during the monitoring periods (Supplementary Fig. 5 Table 3). Under the range of water discharges from 0.1 to 100 m 3 s −1 , river SS carrying capacity exhibited a considerable increase from 2013 to 2016 and a slight decrease after decontamination (2017 to 2018). Contrastingly, the range of water discharges was relatively narrow upstream, and a steady decrease in SS loads has been observed since 2015. Although the above result suggests an increase in SS supply during the decontamination stage, the high SS carrying capacity in 2015 is not consistent with the actual decontamination progress. Since decontaminated regions tend to be bare land, sediment loads are prone to increasing during rainstorms due to soil erosion. Governmental decontamination plans showed that over 50% of agricultural land decontamination was planned to be implemented in 2016, implying that the erosion potential of the decontaminated regions should have been higher in 2016, rather than in 2015.

and Supplementary
The variations in downstream river SS loads over the 6 years (Supplementary Table 4) exhibited a similar trend to peak river SS load in 2015 (126.7 ± 0.3 Gg yr −1 ). This was an approximately 1,776%, 140% and 215% increase relative to that of 2013, 2016 and 2018, respectively. The historical rainfall records (Fig. 4b) show that the rainfall in September 2015 (551 mm) was more than two-fold greater than that during the same period in 2016 (274 mm), implying that the SS peak may be related to strong runoff. Here we estimate the SS loads at 1-year-flood discharge (Q = 95 m 3 s −1 ) using established L-Q curves, which allow for the comparison of dynamic variations in SS loads under the same flood conditions. In Fig. 4b, a significant increasing trend during the decontamination period is shown, with a 237.1% increase in 2016 compared with 2013. After decontamination, the SS loads drastically decreased by approximately 41% from 2016 to 2017, implying changes in sediment yield and transfer patterns due to natural restoration. These results reveal that river SS loads responded closely to land-cover changes during the study period.
To better explore the response of river SS load to land-cover changes, we extracted river monitoring data during each rainfall event and quantitatively linked the river SS to the corresponding soil loss from the decontaminated regions. We found that SS loads during rainstorms were highly correlated with water discharges in both upstream and downstream areas ( Supplementary  Fig. 6). Comparing similar rainfall events, significantly greater SS concentrations are observed in 2015-2016 than in other years (Supplementary Figs. 7 and 8). Considering that the land-cover changes induced by decontamination were more pronounced in the summer season, the regression was performed for SS loads between May and October and soil loss during the corresponding period. A more significant correlation was observed (Fig. 4c) between estimated soil loss by RUSLE and SS load upstream (R 2 = 0.55, P < 0.01, N = 34) than downstream (R 2 = 0.27, P < 0.01, N = 52). Eliminating the effect of rainfall and normalizing by discharge (Fig. 4d) results in a more evident relationship between the erosion potential and SS loads downstream (R 2 = 0.35, P < 0.01, N = 52). Overall, these results demonstrate the connection between river SS dynamics and land-cover changes in the decontaminated regions. The short distance between the upstream catchment and decontaminated regions makes soil erosion a critical driver for upstream river SS transport, whereas the downstream river is dependent on long-distance SS transport, making water discharge an important driver for the downstream catchment. . c, Relationship between the RUSLE-based soil loss (that is, R × K × LS × C × P) from decontaminated regions and river SS load for each rainstorm event (i.e. points) in the downstream (Haramachi) and upstream (Notegami) catchments. d, Relationship between the erosion potential (that is, K × LS × C × P) and SS load (normalized by water discharge) for each rainstorm event (i.e. points) in the downstream (Haramachi) and upstream (Notegami) catchments. The shaded areas represent the 95% confidential interval of the fitting curves (i.e. solid lines).  Cs (30.1 yr), the early decontamination period (16.9 yr) and the surrounding contaminated catchments (mean of 4.92 yr) 20 . Such a sharp decrease in the 137 Cs concentration was also observed at the other three monitoring sites (Supplementary Fig. 9a). Because the 137 Cs concentration in decontaminated soil was considerably lower than that in the contaminated soil 37,42 , these results suggest the contribution of sediment from decontaminated regions to the river system. Moreover, strong negative correlations were observed between measured 137 Cs levels and decontamination progress at all monitoring sites (Fig. 5b and Supplementary Fig. 9b), which further supports our interpretation. The observed 137 Cs concentration increased by approximately 150% in 2018 compared with that at the end of 2016, which may be caused by the weakened sediment supply from decontaminated regions owing to natural restoration and resulting in a relatively increased contribution of sediments from contaminated forest regions 13 .

long-term impact on river SS and 137 Cs discharge
Given that the variation in 137 Cs concentrations reflects a change in sediment source, the deviations between the measured 137 Cs and the natural decrease in 137 Cs derived from surrounding contaminated catchments provide a way to quantitatively estimate the contribution of sediment from decontaminated regions (Fig. 5c). In the early decontamination period, our results showed slight variations in the 137 Cs concentration (Fig. 5a), which could be attributed to the contribution of sediment from upstream regions with different degrees of contamination. During the main decontamination period, the erosion potential in the summer of 2015 was approximately 22% higher than that in 2014 and the heavy rainfall caused the largest flooding event during the study period (26-year flood). This may result in the sediment from decontaminated regions not being the dominant source for downstream. In 2016, decontamination caused an increase by approximately 59% in the erosion potential compared with 2013, with the contribution percentage steadily increasing over this period to a maximum of 75.7% ± 3.2% (value ± 95% uncertainty). After decontamination, the decreased contribution of sediment from decontaminated regions and the increased 137 Cs concentrations can both be attributed to the reduction of soil loss from upstream due to natural restoration.
The 137 Cs discharge from contaminated catchments around the Fukushima region into the Pacific Ocean is another ecological issue of global concern. Our data show that the export flux of particulate 137 Cs from the downstream of the Niida river (that is, Haramachi) peaked in 2015 (1.24 TBq yr −1 , equalling 0.65% 137 Cs loss), which is an approximately 667%, 233% and 429% increase relative to that in 2013, 2016 and 2018, respectively (Fig. 5d). Although such 137 Cs loss is negligible compared with the terrestrial inventory, it is approximately 10 5 times greater than that in the pre-Fukushima stage 43,44 . Accordingly, the dynamic variations in 137 Cs discharge from terrestrial environments into the Pacific Ocean, and its drivers, require more attention in the future.
Here we used SS loads at one-year-flood discharge to normalize the 137 Cs flux and found its peak occurred in 2015 (Fig. 5e). Additionally, the reduction of the normalized 137 Cs flux from 2013 to 2016 (~32%) was similar to natural attenuation in the non-contaminated catchment (~34%) over same period, which may be due to the increased SS load during decontamination offsetting the role of declining 137 Cs concentrations in reducing 137 Cs emission. During the subsequent natural restoration period, the rapid NDVI increase (Fig. 3f) suggested vegetation recovery in decontaminated regions and a decrease in regional erosion potential (Fig. 3g). This resulted in an approximately 24% decrease in sediment yield from the catchment and an approximately 31% decrease in the contribution of sediment from decontaminated regions (Fig. 5c)

Discussion
Our work highlights the great potential of interdisciplinary analyses for understanding river SS variation and quantifying the contribution of sediment from specific regions. Fukushima decontamination practices, like a controllable validation experiment, justified the reliability of using long-term 137 Cs monitoring data for tracing sediment source dynamics due to specific perturbation. Combining the long-term dataset of 137 Cs (or other fallout radionuclides) in SS with remote sensing images would provide additional evidence to determine if the changes in the downstream SS transport pattern are linked to the upstream perturbation.
With these interdisciplinary analyses, we systematically reveal how changes in land use in the decontaminated regions significantly influences downstream river SS and 137 Cs discharge into the ocean. Indeed, the secondary environmental impacts of surface remediation are being increasingly considered in the broader field concerning remediation of regions contaminated with hazardous materials (for example, heavy metals and organic contaminants) 45 . The concept of environmental sustainability-centred green remediation has also been brought up in many scenarios 46,47 . The Fukushima decontamination practice provides evidence showing that mechanical remediation can cause persistently excessive SS load downstream, though it also reduced river 137 Cs concentrations. Since persistently excessive turbidity in rivers affects not only surrounding residents' water use but also trophic level structure in aquatic environments 48 , the unsustainable downstream impacts caused by upstream decontamination should be highly regarded. The vegetation recovery after land development is highly dependent on local conditions 49 , and the soil used for decontamination and local high rainfall amount in Fukushima promoted rapid vegetation recovery 11,13 , which shortened the duration of such unsustainable impacts. Therefore, future upstream contaminated lands that await mechanical remediation need to consider the pre-assessment of local natural restoration conditions or the preparation of appropriate revegetation measures in the catchments' regulatory frameworks, which would minimize the impact of long-term decontamination on downstream sustainability.

Methods
Study region. The Niida River Basin (265 km 2 ) is located about 40 km northwest of the damaged FDNPP. The topography of its upstream is almost mountainous and its soil types are mainly cambisols and andosols, while fluvisols are the dominant soil type in the downstream plain 50 . The monitoring data from the Japan Meteorological Agency show that the average rainfall in the Niida River Basin is greater than 1,300 mm, with more than 75% of the rainfall occurring between May and October. According to the third airborne monitoring survey by the Japanese government, the 137 Cs inventory in the Niida River Basin was over 700 kBq m −2 (ref. 14 ). Because of particularly high contamination in its upstream watershed (over 1,000 kBq m −2 ) 51 , the government-led decontamination was implemented in the upstream basin from 2013 to 2016 (~1% of the area was extended to March 2017).
Land-cover observation. We constructed the vector decontamination maps based on the paper maps from the Ministry of the Environment, Japan. The boundaries of the decontaminated regions with different land-use types were first outlined by creating polygons using Google Earth. Subsequently, the projections of these polygons were imported to ArcMap v.10.3 to quantitatively evaluate their area.

Quantification of land-cover changes in decontaminated regions.
We calculated NDVI within the boundary of the decontaminated regions to quantify the land-cover changes. Through the spectral reflectance dataset in the red (R, nm) and near-infrared (NIR, nm) regions, the NDVI was calculated as 52 : The available satellite images from 2011 to 2018 from Sentinel 2 were downloaded from the United States Geological Survey 53 , while the concurrent MODIS images were derived from the National Aeronautics and Space Administration's Reverb 54 . The wavelength bands and spatiotemporal resolutions of the satellite images used here are summarized in Supplementary  Tables 5 and 6.
To confirm the reliability of the newly generated NDVI variation curve, NDVIs for the same date as the Sentinel 2 images were estimated using the interpolation and compared with the Sentinel 2-based NDVI. The linear regression analysis showed that the fitting R 2 was 0.99 (N = 16, P < 0.01). We also calculated the NDVI in the decontaminated region based on available satellite images of Landsat 5/7/8 (ref. 53 ) and established a daily NDVI variation curve. The linear regression analysis also showed a high R 2 between two daily NDVI curves (R 2 = 0.97, N = 2868, P < 0.01). Therefore, these NDVIs calculated by different satellite images confirmed the reliability of the NDVI variation curve based on ESTARFM.

Estimation of erosion potential in decontaminated regions.
To link the land-cover changes in decontaminated regions with the soil erosion dynamics, we defined an erosion potential (K × LS × C × P) based on the RUSLE.
The soil loss (A, t ha −1 yr −1 ) of a specific region can be estimated as 55 : where R is the precipitation erosivity factor (MJ mm ha −1 h −1 yr −1 ), K represents the soil erodibility factor (t h MJ −1 mm −1 ), L and S are slope length factor (dimensionless) and slope steepness factor (dimensionless), respectively, and C and P are the cover management factor (dimensionless) and support practice factor (dimensionless), respectively. Because these parameters are often set as fixed values, it is difficult to assess the soil loss dynamics during anthropogenic disturbances. To address this problem, we used daily NDVI data to estimate C × P and then considered these dynamic factors in RUSLE. Wakiyama et al. 40 reported a correlation between vegetation cover in Fukushima and the sediment discharges from the standard USLE plot (that is, soil loss, A) that have been normalized by R, K, S and L factors 40,56 . Therefore, this empirical equation reflects the quantitative relationship between vegetation fractions (VF) and C × P.
To quantify daily C × P changes in decontaminated regions, we first converted the interpolated daily NDVI into the VF by a semi-empirical equation 41 : where NDVI s and NDVI ∞ represent the NDVI value for land cover corresponding to no plants and 100% green vegetation cover, respectively. Since these values mainly depend on plant species and soil types, we followed previous methods applied to agricultural land and set NDVI s and NDVI ∞ as 0.05 and 0.88, respectively 41 . Subsequently, the C × P was estimated by the empirical equation derived from uncultivated farmlands and grasslands (R 2 = 0.47, N = 145) 40 : Since the soil type used for decontamination is generally the same, the K factor was set as a constant (0.039; ref. 40 ). For the LS factor, we downloaded a digital elevation model from the Geospatial Information Authority of Japan (spatial resolution: 10 m) to build an LS map using 55 : where Q a is the flow accumulation grid, S g represents the grid slope as a percentage, M is the grid size and y is a parameter depended on slope steepness. We here used the y values recommended by a published study, ref. 55 . The calculated LS-factor map ( Supplementary Fig. 2) showed a relatively consistent LS distribution in space. Based on the ESTARFM-generated satellite images, we compared C × P and erosion potential (K × LS × C × P) and found a significant correlation (R 2 = 0.99, P < 0.01, N = 174). Since these results suggest that LS factors in decontaminated regions have a negligible effect on the erosion potential, the mean LS factor and interpolated NDVI based on the daily variation curve (Fig. 3b) were used to estimate the daily erosion potential.
Monitoring of river discharge and turbidity. The water-level gauges (in situ Rugged TROLL100 Data Logger) and a turbidimeter (ANALITE turbidity NEP9530, McVan Instruments) were installed in each monitoring site to continuously recording the water level and turbidity with a temporal resolution of 10 min. As ocean tides may influence the accuracy of water-level monitoring, the Sakekawa site (M4 in Fig. 2) was excluded from the river monitoring programme.
The recorded water level (H, m) was converted to the water discharge (Q, m 3 h −1 ) based on the annual H-Q curves for each monitoring site. These curves were calibrated using a synchronous monitoring dataset of 10-min-resolution water level and discharge provided by the Fukushima prefecture's official monitoring network 57 . Because of occasional damage to the water-level gauge at the Haramachi site, the available monitoring data with a temporal resolution of 10 min recorded by the Fukushima prefecture's official monitoring network 57 were used to fill the gaps. The percentages of filling data from official monitoring network were all less than 34% except for 2015 (56.6%). Although similar situations occurred in Notegami, we were unable to fill in gaps with other data due to the lack of a concurrent monitoring network.
The hourly SS concentration (C ss , g m −3 ) at each monitoring site was calculated from the measured turbidity (T, mV) using a calibrated curve 20 . As the turbidimeter was susceptible to the moss and debris flowing in the river, the dataset was verified with an automated check by HEC-DSSVue (The U.S. Army Corps of Engineers' Hydrologic Engineering Center Data Storage System) before transforming the data.
The SS load was estimated as the product of the corresponding datasets of discharge and SS concentration, after which we can obtain the annual SS load (L, ton yr −1 ) by taking the sum: We estimated values for gaps including missing and abnormal data through a linear model established by 10-min-resolution monitoring data at the same site. The reliability of the gap-filling strategies used in this study has been documented by Taniguchi et al. 19,20 These procedures vastly enhance the possibility of reconstructing the complete dataset. In this study, only the error in converting from water discharge to SS load was considered in the uncertainty assessment, and all estimates were within 0.5% (95% confidential interval) in this case. To reduce the uncertainty of L-Q fitting, the 10-min monitoring dataset (discharge and SS load) was transformed to a 1-hour dataset.
Considering the river SS is often transported by discharge, we used downstream L-Q curves to estimate river SS loads at 1-year-flood discharge (Q = 95 m 3 s −1 ), which eliminates the influence caused by different annual water discharges. The 1-year-flood discharge was calculated from the daily maximum discharge data from 1 January 2013 to 30 September 2020 at the Haramachi site.
To compare river SS dynamics during rainfall events, we here defined a rainfall event as the increase in water discharge exceeding 1.4 and 1.6 times the baseflow before precipitation for the upstream and downstream catchments, respectively. As a result, a total of 64 and 72 rainfall events from the Notegami and Haramachi sites were identified.
To study the dynamic relationship between soil loss from decontaminated regions and river SS load, we estimated eroded soil amount during each rainstorm using RUSLE. Specifically, the NDVI during a specific rainfall was determined by interpolation. Subsequently, the corresponding C × P can be estimated using equations (3) and (4). With the mean values of the K and LS factors, the erosion potential can then be calculated. Finally, precipitation erosivity factor (Supplementary Table 7) for each rainfall event can be calculated as 58 : where n is the number of years used, m j is the number of precipitation events in each given year j and E and I 30 represent each event's kinetic energy (MJ) and maximum 30 min precipitation intensity (mm h −1 ), respectively, for each event k. The event's erosivity, EI 30 , can be calculated as 58 : where e r denotes the unit rainfall energy (MJ ha −1 mm −1 ) and v r provides the rainfall volume during a set period (r) (mm). For this calculation, the criterion for the identification of a precipitation event is consistent with previous work, that is, the cumulative rainfall of an event is greater than 12.7 mm (ref. 58 ). If another rainfall event occurs within 6 h of the end of a rainfall event, they are counted as one event. Therefore, the unit rainfall energy (e r ) can be derived for each time interval based on rainfall intensity (i r , mm h −1 ) 58 : The calculation's required parameters were derived from the historical precipitation record from the Japan Meteorological Agency 59 . For the Notegami catchment, the precipitation monitoring data were derived from Iitate. For the Haramachi catchment, the precipitation was obtained from three adjacent monitoring sites (that is, Haramachi, Iitate and Tsushima) with the specific weights of 0.143, 0.545 and 0.312, respectively. These weights were determined by the Voronoi diagram method in a Geographic Information System 60 .
River monitoring of particulate 137 Cs. At each monitoring site, the suspended sediment sampler proposed by Phillips et al. 61 was installed at 20-30 cm above the riverbed for the time-integrated sampling of river suspended sediment. The reliability of this sampler has been widely proven in past studies 19,20 . After sampling, the trapped turbid water and SS samples were transferred into a clean polyethylene container and stored until laboratory analysis.
The SS samples were separated from the collected water mixture via natural precipitation and physical filtration, dried at 105 °C for 24 h and subsequently packed into a plastic container. The activities of 137 Cs in the SS samples (C, Bq kg −1 ) were determined via the measurement system, which consists of a high-purity germanium γ-ray spectrometer (GCW2022S, Canberra−Eurisys, Meriden) coupled to an amplifier (PSC822, Canberra, Meriden) and multichannel analyser (DSA1000, Canberra, Meriden). The measurement system was calibrated with the standard soil sample from the International Atomic Energy Agency. Under the 662 keV energy channel, each measurement batch would take approximately 1-24 h to make the analytical precision of the measurements within 10% (95% confidential interval). All measured 137 Cs concentrations were decay-corrected to their sampling date. Moreover, the results obtained in this study were also normalized by their initial average 137 Cs inventory in the catchment (D, Bq m −2 ) to eliminate the effect caused by spatial differences.
As 137 Cs concentration in the sediment sample depends on particle size 19,20 , we conducted a particle size correction for all measured data in Takase, Ukedo and Haramachi to eliminate this effect. The particle size distributions for dried SS samples were analysed using the laser diffraction particle size analyser (SALD-3100, Shimadzu Co., Ltd.). With the parameterized particle size distributions, the particle size correction coefficient (P c ) can be calculated by 19 : where S r and S s represent the reference and collected samples' specific surface areas (m 2 g −1 ). The exponent coefficient, v, is a fitting parameter associated with chemical and mineral compositions. In this study, the same parameters measured in the Abukuma River, the major river in the Fukushima area, were applied for S r (0.202 m 2 g −1 ) and v (0.65). The specific surface area for collected SS samples was estimated by the following equation under a spherical approximation 20 : where ρ is the particle density and d i and p i denote the ratio and diameter of the particle size fraction for particle i. Therefore, the 137 Cs concentration corrected for particle size can be obtained by dividing the measured 137 Cs concentration by P c .
Considering that the decrease in particulate 137 Cs concentration in a catchment was also affected by natural attenuation, there is a need to eliminate this effect from the declining trend of our observed 137 Cs dataset to highlight the impacts of decontamination. The Ukedo and Takase are rivers surrounding the Niida River with similar contaminated situations. Our long-term 137 Cs monitoring data from downstream of these two catchments showed that their 137 Cs decline trends were relatively steady. Although there is a dam reservoir upstream of Ukedo, the 137 Cs concentration observed both upstream and downstream showed a similar declining trend 62 . Therefore, the above evidence suggests that natural attenuation was the dominant factor controlling the 137 Cs decrease in these two rivers.
Here we assumed that the natural attenuation trend of 137 Cs in the surrounding catchments (Ukedo and Takase) with little effect by decontamination was similar to that of the Haramachi catchment. Thus, we fitted their time change curves of 137 Cs concentration (normalized by average 137 Cs inventory of the corresponding catchment) using an exponential model. We then estimated the 137 Cs concentration at the same sampling time as Haramachi in the two catchments by using the fitting models. Finally, we calculated the mean value of the two datasets and recalculated the effective half-life (T eff = ln(2)/λ; λ is the fitting exponential term) of the natural attenuation by the exponential model.
The 137 Cs flux (L Cs , Bq) for each monitoring site was estimated by the product of the SS flux and the 137 Cs concentration in the suspended sediment sample. We then took the sum over that year: According to the law of error propagation, we considered the error from SS load and 137 Cs measurement in the combined uncertainty assessment for 137 Cs fluxes and found their values are all within 1.1% (95% confidential interval).
Using 137Cs as a tracer in estimating SS source contribution. Although 137 Cs has been widely used in tracing sediment source, the spatial variability of the 137 Cs deposition inventory in the Fukushima catchment hinders the estimation of the source contribution from a specific region. However, for the decontaminated catchments, as the 137 Cs concentration in decontaminated soil was much lower than that in contaminated regions (for example, forested regions and the riverbank) 11,13,63 , the fluctuations in the particulate 137 Cs concentration can help to identify the sediment from the decontamination regions. Specifically, we assumed that the particulate 137 Cs concentrations in surrounding contaminated watersheds (that is, having similar land-use composition) follow a similar decline trend driven by natural reasons, while the decontamination-induced land-cover changes cause other sediment sources to mix with the original river SS and thus result in a deviation in observed 137 Cs concentrations from this natural trend. Therefore, the relative contribution (RC) of the specific sediment source can be expressed as: where the C m is the measured 137 Cs concentration and C s and C n represent the 137 Cs concentration in a specific sediment source and the naturally varied 137 Cs concentration at the same time as the measured 137 Cs. For data comparability, all 137 Cs concentrations presented here were corrected by their particle size and 137 Cs inventory. We also excluded the samples with collection weights below 0.5 g from the calculations due to their high uncertainty in P c measurement. In this study, the specific sediment source is the decontaminated soil in the Niida River Basin where the 137 Cs concentration was approximately 53.99 ± 40.90 Bq kg −1 (refs. 37,42 ; mean ± standard deviation, N = 8). The natural decline of the 137 Cs concentration (that is, λ) was established using temporal variation in 137 Cs data originating from the Ukedo and Takase rivers, which were scarcely influenced by decontamination. The first measured 137 Cs data in Haramachi were set as the starting point of its natural decline curve. The total uncertainty for the contribution percentage of SS from the decontaminated regions was calculated by the propagation of error from each part with the uncertainties originating from the measured 137 Cs concentration, 137 Cs concentration in a specific source and the natural 137 Cs concentration. For the uncertainty in the 137 Cs concentration in decontaminated soil (C s ), we set the standard deviation as its error source, while the natural 137 Cs concentrations were calculated by the propagation of the 95% confidential interval of the fitting curves.
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.