Updated trends of water management practice in the Italian rice paddies from remotely sensed imagery

ABSTRACT The rice growing district in northwestern Italy, where paddies were traditionally flooded throughout spring, was interested by a general decrease of standing water presence caused by the adoption of dry seeding crop practices, with consequences for water management and for the ecology of breeding waterbirds. This communication analyses changes in flooding dynamics in the last four years, estimating them from MODIS data and comparing results with previous knowledge of the same study area. Results highlighted an intensification of the phenomenon in the north-western regions (–3.3 ± 0.6% per year in the period 2013–2021) and the almost complete loss of flooded surfaces east to the Ticino river (reaching in 2021 5% of the flooded extension estimated in 2000). Such findings highlight the importance of monitoring this phenomenon – considered by other authors as the biggest anthropogenic change in surface water of all Europe since 2000 – in near real time from remotely sensed data to monitor dynamics and support sustainable management of water usage at the district level.


Introduction
The Italian rice district, located along the Po, Sesia and Ticino rivers, is not only the largest intensive rice cultivation area in Europe (FAOSTAT, 2013) but also the location in which the largest anthropogenic change of surface water in Europe is occurring (Zampieri et al., 2019). In fact, after 1990, the traditional rice cultivation practice, which involves permanent flooding of cultivated surfaces, is being joined by a new technique in which rice paddies are not flooded until the 2-4 leaf stage. This practice -dry seeding with delayed flooding -is proposed for its agronomical and environmental advantages: it has the benefits of reducing the costs of water purchase and farming equipment (Cesari de Maria et al., 2017), it offers an overall reduction of the global warming potential of rice wetland (net balance between decreases of methane emission and increase of nitrogen oxide; Zampieri et al., 2019) and it provides a more efficient use of seeds (Ente Nazionale Risi, 2019a). The adoption of this practice was limited in the past because of the exposition of seeds and early stage plants to nocturnal thermal stress (Sipaseuth et al., 2007); thanks to the warmer temperatures recorded in spring, this problem is less and less a limiting factor, allowing dry seeding becoming increasingly adopted (Ente Nazionale Risi, 2015;Zampieri et al., 2019). In this scenario, monitoring the distribution of floodings and the diffusion of different rice crop cultivation practices is becoming even more crucial for both researchers and irrigation authorities. Remotely sensed data and methodologies are a valuable instrument to monitor water dynamics; in particular, coarse resolution datasets characterised by a high temporal frequency were broadly used for this purpose, in particular, for extremely dynamic and human-related phenomena such as irrigation practices (e.g., Ahamed & Bolten, 2017;Boschetti et al., 2009;Kuenzer et al., 2015;Lopez et al., 2020;Peng et al., 2011).
A recent study (Ranghetti et al., 2018) discussed the agronomic and environmental advantages of the rice dry seeding practice, comparing it with the problems that an excessive spread could cause, and applied a method  to estimate the extension of flooded surfaces from MODIS data (NASA LP DAAC,) to analyse the evolution of water-and dry-seeded rice fields from 2000 to 2016 in the Italian rice district, highlighting a generalised decreasing trend of flooded surfaces (−0.86% per year). The authors underlined the importance to continue monitoring this phenomenon in the subsequent years in view of the critical water reduction already recorded at the end of the analysed time window. Recent evidences of an additional intensification of the phenomenon were reported, 1,2,3 sharpening the importance of doing that; moreover, local authorities warned about the possibility that an uncontrolled increase of dry seeding would lead to a surplus of water request from late June, when water availability is lower and request is higher due to concurrent crop requirements (Cesari de Maria et al., 2017;Romani, 2008), which cannot be satisfied by the reservoir capacity of irrigation networks for the whole rice district in the case of additional increase of dry-seeded rice crops .4,5 In addition, recently, the research community has considered the impact of dry seeding on water saving at the district level: despite the positive contribution already demonstrated by studies based on single-field experiments (Miniotti et al., 2016), new findings highlighted that the reduction of water floodings on a vast area can negatively impact the groundwater level (Mayer et al., 2019).
Moreover, additional consequences of dry seeding practice are the increase of nitrous oxide emissions, which could be three times higher than using permanent floodings (Kritee et al., 2018;Zhou et al., 2020), and the reduction of the breeding habitat of herons and egrets, which uses these standing water surfaces, a surrogate of natural wetland, as the favourite foraging area. Studies have demonstrated that populations showed a decreasing trend in the last few years and that this effect is likely related to the loss of favourable habitat (Fasola et al., 2014;Fasola & Cardarelli, 2015;Imperio et al., 2017;Ranghetti et al., 2018).
Under this framework, this communication aims to update our knowledge about the variation in irrigation dynamics by applying the above-mentioned method (Ranghetti et al., 2018) to the period 2017-2021, in order to verify how irrigation practices for flooded rice during planting are evolving and whether decreasing trends are continuing.

Study area
The study area corresponds to the largest Italian rice cultivation district, located in northwestern Italy (Piedmont and Lombardy regions; see Figure 1 and {extcolor} 2, part a). In the northwestern portion of this area (provinces of Novara and Vercelli), rice is the dominant crop, whilst in the southeastern one (Pavia and Milano), rice fields are more mingled with other crops (corn, soybean, wheat, barley and permanent meadows). Rice is generally seeded from the beginning of April to the end of May (depending on the variety), whilst harvesting is performed in late September-October (Azar et al., 2016).
In order to analyse local differences in terms of water management, the area was split into 7 sub-districts, taking into account the boundaries of the five local irrigation authorities (Consorzi d'irrigazione) and administrative entities (provinces). The boundaries used in Ranghetti et al. (2018) were maintained in this work.

Satellite data
MODIS TERRA Surface Reflectance data (MOD09A1 500 m dataset; Vermote, 2015) were used to update the multitemporal maps of NDFI (Normalized Difference Flood Index; Boschetti et al., 2009) and NDVI (Normalized Difference Vegetation Index; Tucker, 1979) already produced in Ranghetti et al. (2018) for the period 2000-2016; updated time series were produced for a period of 22 years . Although the quality of MOD09A1 images is generally high (being an aggregate product), pixels characterised by low or medium quality (presence of clouds or cirrus, cloud proximity) were filtered out. For this work, a total number of 330 MOD09A1 images were used (15 images

Reference data
Information about farmers' voluntary declarations in the period 2016-2019, collected by the Italian organisation Ente Nazionale Risi (ENR), was used for validation purposes (see Section 2.4). In the provided dataset, each declaration includes the cultivated variety, the size of paddies seeded in dry condition and the total rice size. Declarations were first aggregated at the municipality level and then at the sub-district level in order to be compared with remote estimates.

Data processing and availability
MODIS data were processed following the methodology exposed in Ranghetti et al. (2018) and schematised in Figure 1 in order to obtain the following yearly output products: FF avg , the averaged Flooding Fraction (FF) in the considered time window (March to June) of each year, is a raster output (map with a resolution of 1 × 1 km) computed aggregating all the FF maps computed from MODIS images during the sowing period of each year, used as a proxy of the overall flooding conditions during the sowing period; WS, the proportion of Water-Seeded rice surface for each year, is a tabular product (one time series over each subdistrict and one over the whole district) obtained by weighting the yearly maximum FF by the proportion of arable land cultivated with rice and calibrating this estimate with reference data in the period 2004-2015. The following equation was derived from calibration and used to compute WS: with WS raw being the uncalibrated WS measures obtained from MODIS data.
In order to ensure that Equation (1) can also be used for the period 2017-2021, in this study, WS estimates were validated against reference data provided from ENR as described in Section 2.3 for the period 2016-2019. To assess derived model performance and validity, the following metrics were considered: model EFficiency coefficient (EF), Root Mean Square Error (RMSE), Mean Absolute Error (MAE) and Mean Bias Error (MBE).
These two datasets were publicly released under a Creative Commons 4.0 license (Ranghetti, 2021).

Data analysis
Data analysis conducted in this work was finalised to update the 2000-2016 analysis and test how temporal trends of FF avg and WS observed in Ranghetti et al. (2018) changed in the last four years.
To do so, an aggregation over 5-year time windows (2001-2005, 2006-2010, 2011-2015 and 2016-2020) was first performed. This operation allowed a descriptive comparison of the evolution of these metrics over time, without the influence of yearly fluctuations.
To analyse the significance of these results, FF and WS were each regressed against year for each of the 7 sub-districts and for the whole rice district.
In order to check if the trends changed at a certain date, a breakpoint analysis (Muggeo, 2003) was performed over the whole district and each sub-district. The existence of a breakpoint in each time series was checked by testing regressions 2 and 4 using Davies' tests (Davies, 1987); the position of the breakpoint and the slope of the regression before and after it were thus quantified by building segmented relationships in the cases a significant breakpoint was found (with a starting value of breakpoints set to 2010 -the central year of the analysed period -and 20 bootstrap samples). Figure 1. Schematisation of the processing steps performed to obtain FF raw and WS products. Thumbnails show which products were obtained as the raster image or municipality/subdistrict vector aggregation and which products are composed of daily images (first three datasets) or yearly seasonal aggregations.
The entire process of data collection, processing and analysis was performed using free software. 6

Validation of the WS product
The validation of WS values against reference data, performed as described in Section 2.4, resulted satisfying EF = 0.97, RMSE = 0.04 and MAE = 0.03. These metrics show a high predictive power of the model (high EF) and good estimation capacity with values even better than the ones obtained in Ranghetti et al. (2018) (EF = 0.94, RMSE = 0.06 and MAE = 0.05 using a leave-one-out cross-validation scheme). Moreover, estimates did not result biased (MBE = -0.02, lower than MAE). For these reasons, the parameters of the logistic regressions were not updated with the ENR new dataset and the model was considered robust to estimate the WS variable in the time window object of this work (2017-2021).

Spatiotemporal distribution of FF avg and WS
The spatial distribution of the averaged flooded surface (FF avg ) is shown in Figure 2. Part (b) shows aggregated values over 5-year time windows for the period 2001-2020: in this way, the changes occurred in the last 5 years, specific topic of this work, can be easily compared with the situation of the previous periods. Table 1 quantifies the areas of flooded surfaces during rice seeding (obtained by multiplying values of FF avg by the rice area of each sub-district) aggregated on the same time windows, together with the 5-year averaged proportion of water-seeded rice surfaces (WS).
These results depict a clear decrease of FF avg and WS during the last 5 years, which affects all the sub-districts. Within the north-western portion, the decrease in the 5-year window 2016-2020 was stronger with respect to the variation of the 5-year period 2011-2015, particularly in sub-district B (Novarese) and the eastern portion of C (Grange vercellesi). Since the north-western subdistricts are characterised by greater flooded surfaces, the reduction of water seeding within these regions drove the decline of the overall water-seeded surface, which reached only 44% of the total rice cultivated area (it was 58% in the previous 5-year window). Nevertheless, the decrease of the flooded surface can also be noticed in the south-eastern portion, where a marked reduction had already occurred in the previous 5-year window: within sub-districts F (Sud Milano) and G (Pavese), the decrease has more than halved, from 35 km 2 in 2011-2016 to 14 km 2 in 2016-2020.
Yearly metrics can give additional details to this scenario. Comparing Figure 2 part (c) and Table 1 part (b) with the same results for the period 2000-2016 (Ranghetti et al., 2018, Figure 7 and Table 5), it is evident that surface areas and WS values for years 2017 to 2021 were lower than in any of the records pre-2015. A slight countertrend can be noticed for years 2018 and 2019: this can be justified by the heavy rainfall events that occurred at the middle of April of these two years, which forced farmers to adopt water seeding practice in late-seeded fields (Ente Nazionale Risi, 2018, 2019). Conversely, 2020 and 2021 values show a strong reduction affecting all the sub-districts, with a total flooded surface of 419 km 2 (34% with respect to the surface flooded in 2000). Particularly evident is the almost complete abandonment of flooding practices in sub-districts F and G: in 2021, only 5% of the original flooded surface was measured (it was the 22% in 2016), for a global flooded area of less than 5 km 2 .

Analysis of temporal trends
The decreasing trend observed in Ranghetti et al. (2018) not only continued in the subsequent four years but also with an increasing intensity.
Results of the linear regressions of FF avg over time, summarised in Table 2 part (a), confirmed this evidence, already noticed from Figure 2. The slope of the regression over years 2000-2021 (representing the yearly decreasing trend) was globally quantified as a 1% decrease (it was 0.9% over years 2000-2016); this reduction affected 5 of Table 1. Quantification of the flooded surfaces over years and sub-districts (in km 2 ), together with WS values (%): averaged values over 5-year windows -part (a) -and update for years 2017-2021 -part (b  the 7 sub-districts. The variation is higher in the northwestern sub-districts, characterised by wider water surfaces; the only 2 sub-districts in which the decrease did not accelerate are the two south-eastern ones, in which the reduction of the flooded surfaces started earlier and already reached low FF avg values. An even more serious intensification of the changes in the slope of the regressions was evidenced by results of the linear regressions of the percentage of rice surfaces seeded in water (WS) over time (Equations 2 and 4). As Figure 3 and Table 2 part (b) underline, the yearly overall slope was -1.5%, 60% more intense than as previously quantified (-0.9%). This decreasing condition, affecting all the 7 sub-districts, was mainly due to the intensification of changes that occurred in the north-western ones, in which the relative reduction of the regression slopes exceeded 100% in 2 cases: it is the case of sub-districts B and D (Oltrepo casalese), already evidenced in Section 3.2.
Davies' tests confirmed the significance of the existence of a breakpoint in the overall time series (p-value < 0.01), estimated for the year 2013 ± 1. As shown in Table   Figure 3. Temporal distribution of WS values along year and sub-districts. Points are predicted values; blue lines represent the linear regressions over the periods 2000-2021 (solid) and 2000-2016 (dashed); red lines are the segmented regressions (with their 95% confidence intervals -grey bars), whose breakpoint positions are indicated by vertical segments (whiskers represents 95% confidence interval). Figure 3, a significant breakpoint was also found in 5 of the 7 considered sub-districts (the north-western and central ones), corresponding to year 2013 (B, D and E) or 2008 (A and C), whilst in southeastern subdisticts F and G, no breakpoints were found. These findings are coherent with our knowledge about the introduction of the practice of water seeding, which had been introduced before 2000 in the rice cropping area to the east of the Ticino river (Romani, 2008), whilst it became an even more common practice only since 2010 in the western part (Ente Nazionale Risi, 2014). Table 3 highlights another interesting finding. Ranghetti et al. (2018) discussed the differences in the slopes of the WS regressions signalling the importance to monitor if the WS proportion reduction would have progressed in the northwestern area with the same intensity detected for the southeastern one. The breakpoint analysis allowed quantifying the slope of the regressions after the diffusion of the water seeding, showing the intensity of the yearly decrease being between -2.7% (sub-districts A -Baragge -and C) and -3.5 to -3.7% (B and E -Lomellina). Sub-district D showed an anomalous decrease (-5.1%): this is a small region, in which changes can occur suddenly with respect to other ones (moreover, Figure 3 shows that the anomalous values of year 2013 and 2014 might have concurred to raise the breakpoint y value, causing a decrease of the slope). These results demonstrate that, despite global slopes measured over north-western sub-districts being lower than over south-eastern ones, the decrease of waterseeded rice surfaces after the adoption of this practice is even faster here.

Conclusions
This communication updated the estimations of the Flooding Fraction (FF avg ) and the proportion of Water-Seeded rice surfaces (WS) from MODIS data in the Italian rice district in the period 2000-2021, an area considered with the biggest anthropogenic change of surface water of all Europe since 2000. In comparison to the previous analysis (Ranghetti et al., 2018), which had focused over the period 2000-2016, the new results underlined how the reduction of the rice surface flooded during seeding continued in the last four years amongst all the considered sub-districts at accelerated pace. The overall yearly decrease was estimated to be -1.02 ± 0.13% for FF avg and -1.5 ± 0.2% for WS, respectively, 18% and 60% lower than over years 2000-2016, causing the loss of 812 km 2 (66%) of flooded surfaces with respect to 2000 estimations.
In sub-districts located east to the Ticino river, the reduction continued with a constant trend with respect to previous years, reaching in 2020 and 2021 the almost complete loss of flooded surfaces (less than 5 km 2 of flooded surfaces, representing 5% compared to 2000). Contextually, in the north-western regions, where water availability is higher, the reduction after the adoption of dry-seeding practices (that was estimated using a breakpoint analysis to be occurred in the time window 2008-2013) resulted particularly fast (-3.3 ± 0.6% per year): supposing that this trend will continue with this intensity in the coming years, water availability in summer through irrigation networks will not be sufficient to satisfy conflicting crops demand caused by a wider diffusion of dry-seeded rice surfaces.
Findings of this work, together with the necessity not to reduce water-seeded rice surfaces in a critical extent both for reservoir capacity and ecological reasons, underline the importance of a real-time monitoring of flooding conditions, which this method can ensure until MODIS data are provided. For this reason, the authors of this work are planning to disseminate yearly updated data (maps, time series and trends) using an interactive platform, so granting to researchers and irrigation authorities the availability of an updated overlook of flooding condition in near real time, plausibly by the end of June of each year. Furthermore, this method could be adapted to different temperate rice cultivation areas for which analysis of trend in water management is an important topic, e.g., in United States, rice is directly seeded adopting dryand water-seeding methods, whose allocation, according to literature, is 67% and 33%, respectively (Kumar & Ladha, 2011). However, these statistics are based on data collected in the 1980s (Hill et al., 1991); in the light of the above, information about the spatial distribution and the trend of diffusion of water-seeding methods could be useful to update our knowledge and to provide early warnings in the case of potential variations, which could cause criticalities for water management.