Assessment of Water Management Changes in the Italian Rice Paddies from 2000 to 2016 Using Satellite Data : A Contribution to Agro-Ecological Studies

The intensive rice cultivation area in northwestern Italy hosts the largest surface of rice paddies in Europe, and it is valued as a substantial habitat for aquatic biodiversity, with the paddies acting as a surrogate for the lost natural wetlands. The extent of submerged paddies strictly depends on crop management practices: in this framework, the recent diffusion of rice seeding in dry conditions has led to a reduction of flooded surfaces during spring and could have contributed to the observed decline of the populations of some waterbird species that exploit rice fields as foraging habitat. In order to test the existence and magnitude of a decreasing trend in the extent of submerged rice paddies during the rice-sowing period, MODIS remotely-sensed data were used to estimate the extent of the average flooded surface and the proportion of flooded rice fields in the years 2000–2016 during the nesting period of waterbirds. A general reduction of flooded rice fields during the rice-sowing season was observed, averaging −0.86± 0.20% per year (p-value < 0.01). Overall, the loss in submerged surface area during the sowing season reached 44% of the original extent in 2016, with a peak of 78% in the sub-districts to the east of the Ticino River. Results highlight the usefulness of remote sensing data and techniques to map and monitor water dynamics within rice cropping systems. These techniques could be of key importance to analyze the effects at the regional scale of the recent increase of dry-seeded rice cultivations on watershed recharge and water runoff and to interpret the decline of breeding waterbirds via a loss of foraging habitat.


Introduction
Habitat loss is one of the most important threats to biodiversity and is regarded as the major cause of species extinction [1].This is especially true for freshwater environments, which host 6% of known animal species despite covering only 0.8% of the Earth's surface [2] and which suffered the most pronounced decline in biodiversity in recent years [3].In Mediterranean Europe, 80-90% of the natural wetlands have disappeared during the last two centuries mainly due to land reclamation [4], making artificial water surfaces-such as paddy rice fields-a valuable surrogate for the lost natural wetlands as potential habitat for aquatic wildlife communities [5], and particularly for waterbirds [6][7][8][9][10][11][12][13].

The Importance of Monitoring Standing Water Dynamics
During the last decades of the 20th Century, rice fields provided the main foraging habitat for several species of waterbirds that nest in the study area, the most abundant of which were the herons and egrets of the family Ardeidae.
The breeding populations of these herons and egrets have shown an increasing trend from the 1970s to the end of the 20th Century, due to favorable climatic and environmental conditions and to a diminished human-related mortality [21].However, recent studies [22,23] showed that the number of nests of the two most abundant species, the grey heron (Ardea cinerea) and the little egret (Egretta garzetta), decreased in 2016 to about half of that of 2000.Their number continued instead to increase in adjacent natural wetlands, or in the small streams of the uplands above 200 m a.s.l.(Figure 1).The total surface of rice paddies in northwestern Italy remained fairly stable (2048 km 2 in 2000, 2139 km 2 in 2010, 1989 km 2 in 2015; data from [15] and Ente Nazionale Risi, unpublished data; see Section 2.3).Moreover, other possible negative impacts on heron populations such as human-induced mortality and egg contamination by environmental pollutants had already declined sharply before 2000 [21].The most probable cause of such a population decrease is therefore the above-described drastic change in rice cultivation practices happening after 1990, which reduced the suitability of rice paddies as foraging habitat for waterbirds.Total number of nests of the two most abundant breeding waterbirds, the grey heron and the little egret, in north-western Italy during the time window analysed in this study, compared to their predominant foraging habitat [22,23].
In this context, the possibility to map how the adoption of the two main seeding practices changed during the last years could be useful in the framework of studies concerning the population dynamics of waterbird populations.
Geospatial information concerning the extent and changes in rice seeding practices is also important for agronomic reasons and in relation to water management.Rice seeding under flooding, while requiring a large amount of water availability (15,000-20,000 m 3 /ha [24]), allows for a better soil drainage, since aquifers rise during flooding.On the contrary, dry seeding requires instead less water and leads to lower evapotranspiration losses [25], but the water loss by percolation is higher, Figure 1.Total number of nests of the two most abundant breeding waterbirds, the grey heron and the little egret, in northwestern Italy during the time window analyzed in this study, compared to their predominant foraging habitat [22,23].
In this context, the possibility to map how the adoption of the two main seeding practices changed during the last few years could be useful in the framework of studies concerning the population dynamics of waterbird populations.
Geospatial information concerning the extent and changes in rice seeding practices is also important for agronomic reasons and in relation to water management.Rice seeding under flooding, while requiring a large amount of water availability (15,000-20,000 m 3 /ha [24]), allows for a better soil drainage, since aquifers rise during flooding.On the contrary, dry seeding requires instead less water and leads to lower evapotranspiration losses [25], but the water loss by percolation is higher, especially during the first flooding events (beginning of summer), when water availability is generally lower [16,17].Moreover, a delayed maximum water requirement (at the beginning of June instead of at the end of April) can cause competition with other crops, with consequent problems for water management at the district scale [26].The possibility to monitor the reductions in both extent and timing of flooding would be crucial to monitor areas in which water availability could be problematic and eventually to address management actions by public stakeholders to prevent water stress problems.

Remote Sensing Contribution for Flooding Detection
Remote sensing data and techniques are widely used to map and monitor the extent of flooded surfaces [27].Most of the studies dealing with the detection of inundated areas use medium-resolution optical imagery such as those acquired by the Landsat series sensors [28], which is characterized by a fairly high spatial resolution (30 m for sensors working since 1982).This is sufficient for mapping water bodies with few or seasonal variations in their extent like rivers [29][30][31], coast lines [32][33][34], lakes [35], small water bodies [36] or natural stable wetlands [37][38][39][40].However, their 16-day revisiting time limits their usefulness for detecting rapid flooding events or fast water dynamics: this is the case for example of disastrous inundations or agronomic flooding, the persistence of which can be restricted to only a few days [41].The recent availability of Sentinel-2 data [42] allows detecting water bodies at fine spatial resolution (10 m) [43], including flooded rice paddies [44,45] with a shorter revisiting time of five days (thanks to the joint use of Sentinel-2A and Sentinel-2B satellites).However, due to their novelty, these data do not allow performing an analysis of recent multiannual trends of the extent of agronomic flooding.
These considerations make the use of coarse-resolution datasets preferable for the analysis of temporal changes in water management dynamics, e.g., [46][47][48][49][50].In particular, MODIS data [51] appear to be particularly suitable for these purposes, due to their fast revisiting time (daily) and the relatively long period of data availability (since 2000).In a recent study, Ranghetti et al. [52] demonstrated that their spatial resolution, whilst not allowing recognizing flooding conditions in single rice paddies, is sufficient for investigating the spatial and temporal variations of flooded surfaces at the regional scale, at least in the first part of the rice growing period.In particular, they showed that multitemporal MODIS-derived NDFI (Normalized Difference Flood Index [53]) images can be used to derive 1 × 1 km resolution maps of Flooding Fraction (FF: the percentage of flooded surface over the pixel area) within rice paddies.These maps are characterized by satisfactory accuracy when compared against reference ground information (R 2 = 0.73, EF = 0.57 at 1 × 1 km resolution), which makes them useful for multitemporal regional mapping purposes.

Aims of the Work
This work aims to investigate the variation in irrigation dynamics occurring in the last 17 years in the rice paddies of northwestern Italy for an area of about 210,000 ha, by applying the aforementioned method of Ranghetti et al. [52] to generate seasonal FF maps.In particular, the objective is first to verify the presence of a temporal trend in the seasonal flooding fraction at rice district and sub-district levels and then to quantify the proportion of water-seeded rice area and its variations across years.Such data are fundamental to understand hydrological dynamics and their impact on fresh water ecosystems in the study area.In particular, wide area maps of standing water presence are useful to assess the potential effects of variations in water management practices at the regional scale on water table recharge and channels'/rivers' water runoff.Moreover, these data, together with information on other aspects of the foraging ecology of waterbirds within the study area [54] could help shed light on the factors determining the recent decrease of their breeding populations.

Study Area
The study area includes the whole main rice district of northwestern Italy, located between Piedmont and Lombardy along the rivers Sesia, Ticino and Po (Figure 2).
The area was delimited using official municipality boundaries, in order to provide results that can be compared with official statistics.Rice is by far the dominant crop in the northwestern portion of this area (provinces of Novara and Vercelli), while in the southeastern one (Pavia and Milano), rice fields are more mingled with other crops including both summer crops (e.g., corn, soybean), winter crops (wheat and barley) and permanent meadows used for forage production.The area derives the water required for rice paddy flooding from five main irrigation networks (for details, see Figure 2 and Table 1).To highlight local differences in terms of water management and their variations, municipalities were therefore aggregated into seven sub-districts, based on the irrigation network to which they belonged (Table 1).Due to their large dimensions, the Est Sesia and Est Ticino Villoresi networks were further divided according to administrative boundaries (provinces).It is worth noting that small areas belonging to other irrigation networks could be included in each sub-district, since the municipalities' boundaries used to delimit the sub-districts do not always exactly coincide with boundaries of the irrigation network.However, this simplification does not affect the overall approach devoted to regional analysis.Land cover data were used both to mask pixels corresponding to non-arable land and to compute the rice-cultivated area of each sub-district (see Section 2.4.2).In particular, detailed information about the extent of arable land was obtained from two high spatial resolution land use maps: DUSAF 4 (Destinazione d'Uso dei Suoli Agricoli e Forestali, version 4) [55] for Lombardy and LCP (Land Cover Piemonte) [56] for Piedmont (characterized by a spatial resolution of 20 m and 10 m, respectively).
The extent of rice paddies within each sub-district was computed from the Corine Land Cover (CLC) maps [57], Version 18.5.Four CLC products are available: CLC1990 (data collected between 1986 and 1998), CLC2000 (1999-2001), CLC2006 (2005-2007) and CLC2012 (2011)(2012).Since rice extent in the different districts shows very limited variations among the four products (see values in Table 1), the averaged values of the four products were considered representative and valid for the whole temporal range of 2000-2016.

MODIS Satellite Data
MODIS TERRA Surface Reflectance data (MOD09A1 500 m dataset [58]) were used to produce multitemporal maps of NDFI (Normalized Difference Flood Index [53]) and NDVI (Normalized Difference Vegetation Index [59]) over the study area for the 2001-2016 period.This product provides radiometrically-calibrated and atmospherically-corrected data composited using over 8-day periods.It is characterized by a spatial resolution of 500 m, a spectral resolution of seven bands including RED (red band, Band 1, 620-670 nm), NIR (near-infrared band, band 2, 841-846 nm) and SWIR-2 (short-wavelength infrared band, Band 7, 2105-2155 nm), and it is available from February 2000.Information about the day of acquisition (DOY) and the acquisition quality (QA) of each pixel is also available at the same resolution.The MOD09A1 compositing algorithm selects for each pixel the value acquired in a 8-day period characterized by the lower atmospheric disturbance (see the MOD09 user guide [58] for further details), making the presence of cloud-covered pixels much less frequent than in daily data (e.g., MOD09GA dataset).For this reason, MOD09A1 time series are less affected by cloud contamination with respect to daily imagery: this is useful for the analysis of aggregated values within rice sub-districts (see Section 2.4.1), to avoid the possibility that the position and timing of clouds unbalance the aggregation of data.
Images acquired between March and June of each analyzed year (15 composite images for each year, for a total of 255 images) were downloaded for the study area and used to compute the NDVI and NDFI indices using the MODIStspR package [60,61].NDVI [59] and NDFI [53] are computed as follows: where RED, NIR and SWIR are respectively the first, second and seventh MOD09A1 bands.The MODIS Surface Reflectance Data State QA Descriptions layer was used to exclude pixels covered by clouds and cirrus, or with a low data quality flag.In addition, only pixels including at least 50% of arable land were considered, in order to limit the analysis to agricultural areas.

Official Data Concerning Rice Cultivation Practices
The Italian organization Ente Nazionale Risi (ENR) collects farmers' declarations regarding the surface cultivated with the different rice varieties (mandatory information) and the corresponding sowing methods (voluntary information).These data are then processed by ENR to compute aggregated values at the municipality level.Information for the 2004-2015 period was obtained from ENR, compared with MODIS-derived estimates, and then used to calibrate a model for the estimation of the proportion of dry-seeded rice fields starting from satellite data (see Section 2.4.2).

Data Processing and Analysis
Data processing and analysis were performed in three main steps: 1. creation of seasonal maps of average (FF avg ) and maximum (FF max ) Flooding Fraction (FF) amount and timing; 2. computation of the proportion of Water-Seeded rice fields (WS); 3. analysis of temporal trends of flooding extent and timing.
Details about these steps are given in Sections 2.4.1-2.4.3.The retrieved variables are useful indicators for assessing different aspects related to water use and standing water presence in the study area.In particular, FF avg is a proxy of the overall flooding conditions during the sowing period, with higher values indicating a longer persistence of standing water in the fields.FF max is used instead as an indicator of the diffusion of water seeding.Since dry-seeded rice is usually flooded only after rice plants are already grown (third to fourth leaf development stage), when optical remote sensing data are no longer able to detect standing water, the identification of a single flooding event during the sowing period can be assumed as a reasonable indicator of water seeding management.Pixels with higher FF max values indicate therefore areas with a larger presence of water-seeded rice fields.Finally, the WS indicator was computed to assess the proportion of water-seeded rice area with respect to the whole rice-cultivated area and to assess its changes in the considered period.

Creation of Seasonal FF Maps for the Sowing Period
For each available MODIS image (i.e., for each date), FF maps were computed from the NDFI spectral index following Ranghetti et al. [52].According to this method, FF estimates can be considered reliable only when NDVI < 0.4 (i.e., during the first stages of rice growth; about up to the start of the tillering phase) because later in the season, the presence of standing vegetation can mask the presence of water.Moreover, aggregation at 1-km resolution helps to improve the correspondence between FF estimates and ground observations.Therefore, FF predictions for pixels with an associated NDVI > 0.4 were removed from the analysis, and the NDVI-screened FF maps were resampled to 1-km resolution (using average values for FF and NDVI and median values for DOY).For additional details about the accuracy assessment of FF maps, see Ranghetti et al. [52].
The time window considered for this computation was the period between mid-April and the end of May (corresponding to MODIS DOYs of composite between 105 and 145).This was selected based on expert knowledge about typical crop calendars and agro-practices in the area [16], as well as on results of previous studies at regional scale on rice phenology [47,76,77].A preliminary analysis of MODIS data was also conducted to confirm this choice and to test screening criteria allowing unreliable observations to be discarded in the computation of FF avg (see Section 3.1).
Maps of flooding conditions within the sowing season of each year were successively created starting from the multitemporal FF maps.Two different maps with the same spatial resolution (1 × 1 km) were then created for each year: where i is indexing the n MODIS composite images available in the sowing period (composite DOY from 105-145), in year y.
Both FF avg and FF max maps were finally aggregated at sub-district and district level, by averaging FF values respectively within sub-districts and in the whole rice district.
Finally, maps with the DOY (Day Of the Year) of maximum seasonal FF (DOY(FF max )) were generated to check if the seasonal dynamics of FF had been changing since 2000 (see Section 2.4.3).

Computation of the Proportion of Water-Seeded Rice Fields
The proportion of Water-Seeded rice-cultivated area for each year (WS) was computed at district and sub-district levels, in order to analyze how this cropping practice has been changing since 2000.
Raw WS values for each year and sub-district were first computed as: where A arable (s) is the total arable area of the considered sub-district s (or of the whole district), while A rice (s) is the rice area in sub-district s, computed as described in Section 2.2.1.WS raw (y, s) values were then fitted against official data provided by Ente Nazionale Risi (see Section 2.3) using a logistic regression analysis (the regression equation is expressed with the modified Wilkinson-Rogers notation for linear models [78], here and in the following): The official proportion of water-seeded rice fields (WS ENR ; see Section 2.3) was obtained from official statistics as the ratio between the total extent of parcels declared as water-seeded within each sub-district, and the corresponding total rice area (A rice ).
The logit transformation was chosen to constrain WS = 0 to 0 and WS = 1 to 1, correcting only for underestimations or overestimations of intermediate values.The accuracy of the regression model logistic was tested performing a Leave-One-Out Cross-Validation (LOOCV) analysis.

Analysis of Temporal Trends of DOY(FF max ), FF avg and WS
The presence of temporal trends of DOY(FF max ), FF avg and WS values were finally tested.In order to avoid problems of spatial autocorrelation and to depict easily the changes occurring in the analyzed period, the trend analysis was performed using the values averaged at sub-district and district values.
At the rice district scale, univariate linear regressions were used: while, at the sub-district level, the effects of sub-district and of its interaction with year were also assessed by performing a multivariate ANOVA analysis: Finally, univariate regressions were performed to identify the multi-annual trends (magnitude and significance of the slope coefficient) for each sub-district considered by fitting 7 separate models: where s is the sub-district considered in each regression.
To evaluate the existence of statistically-significant trends in the fitted regressions, standard statistical metrics were computed.The value and the standard error of the temporal regressor (year) were used to evaluate the magnitude of the trend in the univariate regressions.The significance (p-value) of the t-value of this regressor (estimate divided by standard error) was considered in order to verify the statistical significance of the trend, using p = 0.05 as the threshold value to identify significant trends.The effect of the regressors in the multivariate analysis is analyzed with a standard Analysis of Variance (ANOVA Type II), using the F-value of the Fisher test to check which regressors explain more variance and the related p-value to verify the statistical significance of their effect.Finally, a qualitative analysis of the variation of the trend slopes was performed by graphically comparing the fitted linear regressions (Equations ( 13) and ( 14)) with smoothed tendencies; for this purpose, LOESS (LOcally-weighted scatterplot smoothing regrESSion) curves were fitted using a span value of 0.75.

Analysis of the Seasonal Distribution of Flooding Conditions
The seasonal analysis of flooding fraction and NDVI variations allowed detecting the time window during which flooding events are mainly concentrated, and the period in which their detection can be considered reliable (NDVI < 0.4) and related to the sowing phase.
Figure 3 shows the temporal variations of the distribution of NDVI and FF values obtained from MODIS in the analyzed period.Water presence was low in March and at the beginning of April (first three composite boxes), confirming that flooding generally starts after the 10th of April.Moreover, NDVI distributions highlight that, after the end of May, the number of unreliable pixels (NDVI > 0.4, following results reported in Ranghetti et al. [52]) exceeds that of the reliable ones.This indicates the presence of already grown rice plants, making the FF estimates highly uncertain.
Therefore, analyses of FF and WS temporal trends were based only on images acquired within composite DOYs 105 and 145 (see Section 2.4.1).It is however important to remind that the aim of this study is to assess variations in flooding conditions during rice-sowing, when the crop is not yet well developed and NDVI is generally low.confirms a global change of 0.75 ± 0.31 days per year (t = 2.39, p-value = 0.03).As shown in Figure 5, the flooding period in north-western sub-districts is generally slightly anticipated with respect to the south-eastern ones (higher DOY value, although this is not evident in all years), consistently with the fact that north-western sub-districts, being nearer to rivers, can be the first to receive water when it becomes available in spring.This is confirmed by the ANOVA (Table 2) of the multivariate regression (equation 10), which points out that most of the variance is explained by the "year" variable (F = 28.29) and only slightly by the sub-district (F = 2.61).No difference in the trend between sub-districts resulted from the analysis (interaction between year and sub-district: F = 0.29, p-value = 0.94), even if the south-east ones revealed a slightly bigger change (Figure 5, panel of the standard deviation).

Analysis of the spatial distribution of flooding conditions
The spatial distribution of overall flooding conditions during the sowing phase (FF avg ) was not uniform, with higher values detected in the north-western parts of the study area than in the south-eastern one.
This difference is also highlighted by Figure 7, where two groups of sub-districts can be easily identified.While in the entire rice district FF avg and FF max averages are 0.22 and 0.35 respectively, The analysis of the temporal changes of maximum FF occurrence (DOY(FF max ); see Section 2.4.3)highlighted a general delay throughout the whole rice district.The DOY of maximum flooding extent changed from about 120-130 (first ten days of May) during 2000-2009 to about 130-140 (mid-May) in the last seven years (Figure 4); the univariate linear regression between DOY(FF max ) and year (Equation ( 7)) confirms a global change of 0.75 ± 0.31 days per year (t = 2.39, p-value = 0.03).As shown in Figure 5, the flooding period in the northwestern sub-districts is generally slightly anticipated with respect to the southeastern ones (higher DOY value, although this is not evident in all years), consistent with the fact that northwestern sub-districts, being nearer to rivers, can be the first to receive water when it becomes available in spring.This is confirmed by the ANOVA (Table 2) of the multivariate regression (Equation ( 10)), which points out that most of the variance is explained by the "year" variable (F = 28.29) and only slightly by the sub-district (F = 2.61).No difference in the trend between sub-districts resulted from the analysis (interaction between year and sub-district: F = 0.29, p-value = 0.94), even if the southeast ones revealed a slightly bigger change (Figure 5, panel of the standard deviation).

276
This difference is also highlighted by Figure 7, where two groups of sub-districts can be easily

Analysis of the Spatial Distribution of Flooding Conditions
The spatial distribution of overall flooding conditions during the sowing phase (FF avg ) was not uniform, with higher values detected in the northwestern parts of the study area than in the southeastern one.
This difference is also highlighted by Figure 6, where two groups of sub-districts can be easily identified.While in the entire rice district FF avg and FF max , the averages are 0.22 and 0.35, respectively, the values between Sub-districts A-B-C (northwestern part, FF avg = 0.34 and FF max = 0.53) and D-E-F-G (southeastern part, FF avg = 0.10 and FF max = 0.18) differ considerably.This is due to various causes, among which: (i) the proportion of rice-cultivated area, which is considerably higher in the northwestern parts, leading to higher FF avg (see Table 1); (ii) the greater and earlier adoption of dry seeding in the eastern part, leading to lower FF max values.Considering the altitudinal gradient of this area, which decreases from northwest to southeast, it is also evident how the averaged flooding fraction decreases from the higher, northwestern sub-districts towards the lower, southeastern ones, in accordance with our knowledge about the geographical distribution of water availability in the study area.Northwestern sub-districts, being nearer to rivers, have higher water availability with respect to the southeastern ones.This probably led to a greater spread of dry seeding in the southeastern sub-districts, as already pointed out in Ranghetti et al. [52].The map of the coefficient of variation of seasonal FF avg values (Figure 7, last panel) also highlights that the interannual variation of standing water presence is higher in the eastern sub-districts, where water is less available.
Version February 16, 2018 submitted to Remote Sens.
12 of 23   This is due to various causes, among which: i) the proportion of rice-cultivated area, which is 281 considerably higher in the north-western parts, leading to higher FF avg (see The temporal analysis of the seasonal average flooding fraction (FF avg ) allows identifying and 301 quantifying a decreasing trend in overall flooding conditions during the sowing period (both globally 302 and within each sub-district), and also identifying differences between sub-districts.  Figure 7 finally highlights that Sub-district E (Lomellina) exhibits differences between its northwestern (higher flooding fractions), northeastern (lower values) and southern (intermediate values) parts.This is due to differences in crop fractional cover (rice is dominant in the northwestern and southern parts, while in the northeastern one, other spring crops prevail) and rice practices (water availability is higher in the northwestern part, particularly in the area between Vigevano and Trecate), which make this sub-district fairly heterogeneous.Nevertheless, it was considered as a unique sub-district to simplify interpretation of the results aggregated on administrative boundaries.

Analysis of Temporal Trends in Flooding Conditions in 2000-2016
The temporal analysis of the seasonal average Flooding Fraction (FF avg ) allows identifying and quantifying a decreasing trend in overall flooding conditions during the sowing period (both globally and within each sub-district), as well as identifying differences between sub-districts.
When the whole rice district is considered, the univariate linear regression between FF avg and year (Equation ( 8)) shows a decrease of −0.86 ± 0.20% per year (t = −4.33,p-value = 0.00059).This metric confirms that flooding practice, hence water diffusion, has been steadily decreasing within the study area.
Results of the ANOVA analysis, reported in Table 3, allow ascertaining the effect and significance of the different factors (years and sub-districts) and of their interaction.Although most of the variance is explained by the sub-districts (an expected result considering the differences in rice cover and irrigation regimes among sub-districts; see Figure 6), the effect of year is highly significant.This confirms that the decrease in flooding occurrence highlighted by the previous univariate regression analysis is a common pattern for all sub-districts and was not conveyed by the effect of averaging values from different sub-districts.The significant sub-district × year interaction, even if low, also suggests that the rate of FF avg decrease varied between sub-districts.This is however influenced by the fact that southeastern sub-districts are characterized by a generally lower FF avg : the slope of the trend line can therefore be expected to be lower even if a similar relative interannual decrease is observed.
Results of the regression analysis conducted at the sub-district level (Figure 8 and Table 4) confirm that sub-districts can be split into two groups according to their mean flooding fraction, as already pointed out in the descriptive analysis.In the first group (Districts A, B and C, i.e., the northwestern ones) FF avg were considerably greater than in the second group (Districts D, E, F and G, i.e., the southeastern part), in accordance with the large amount of variance explained by the sub-district predictor in the multivariate regression (12).Table 3 also highlights the presence of a statistically-significant decreasing trend for all of the sub-districts.In summary, even if overall flooding conditions are spatially heterogeneous, their decreasing trend is widespread over the whole rice cultivation area, with a minimum decrease of −0.21 ± 0.09% per year in Sub-district D and a maximum decrease of −1.46 ± 0.28% per year in Sub-district C. When the whole rice district is considered, the univariate linear regression between FF avg and year (equation 8) shows a decrease of −0.86 ± 0.20% per year (t = −4.33,p-value = 0.00059).This metric confirms that flooding practice, hence water diffusion, has been steadily decreasing within the study area.
Results of the ANOVA analysis, reported in Table 3, allow ascertaining the effect and significance of the different factors (years and sub-districts) and of their interaction.Although most of the variance is explained by the sub-districts (an expected result considered the differences in rice cover and irrigation regimes among sub-districts -see Figure 7), the effect of year is highly significant.This confirms that the decrease in flooding occurrence highlighted by the previous univariate regression analysis is a common pattern for all sub-districts, and was not conveyed by the effect of averaging values from different sub-districts.
The significant sub-district × year interaction, even if low, also suggests that the rate of FF avg decrease varied between sub-districts.This is however influenced by the fact that south-eastern sub-districts are characterised by a generally lower FF avg : the slope of the trend line can therefore be expected to be lower even if a similar relative interannual decrease is observed.
Results of the regression analysis conducted at sub-district level (Figure 8 and Table 4) confirms that sub-districts can be split in two groups according to their mean Flooding Fraction, as already  Values of FF max multiplied by the rice area of each sub-district can be used to quantify the surface that was flooded at least once per year during sowing, a measure that can be considered a proxy of the water-seeded area (Table 5).Estimated flooded surface has reduced from 2000-2016 in all sub-districts, with a global reduction of 44% and a reduction within sub-districts varying from 24% in Sub-district A to 80% in Sub-district G.This finding is particularly interesting from an ecological point of view: the higher relative flooding decrease, i.e., the loss of foraging habitats for herons and egrets, appears in fact to be located within the eastern sub-districts where also landscape fragmentation is more pronounced, with a potential combined negative effect for waterbird populations.
Table 5.Quantification of the surface (in km 2 ) of flooded surfaces over years and sub-districts.Results of the ANOVA analysis on the proportion of water-seeded rice fields (Table 6) shows that, similarly to what observed for FF avg , sub-district is the most important factor in explaining WS variance, but year also has a statistically significant effect.However, differently from what observed for FF avg , the interaction between district and year is not statistically significant (F = 1.19, p-value = 0.32).
The univariate linear regression over years showed a significant WS decreasing trend of −0.91 ± 0.28% per year (t = −3.27,p-value = 0.0051).This decrease is similar to that of FF avg (−0.86 ± 0.20% per year), and corresponds to a reduction of −14.6 ± 4.5% (−376.8± 1.6 km 2 ) throughout the whole study period.The analysis within sub-districts (see Table 7 and Figure 10) shows that, also in the case of WS, a general decreasing trend is observed in all sub-districts, although some are not statistically significant (i.e.sub-districts A, B and D).
However, comparing slope coefficients reported in Tables 4 and 7, it can be observed that, while in the case of FF avg the greater decrease occurred in the north-western sub-districts, in the case of WS the situation is opposite.This apparent contradiction is justified by the fact that the slope of FF avg reduction depends also on the different rice fractional cover.A similar relative reduction of FF avg leads in fact to a larger slope of the regression for sub-districts characterised by a higher FF avg at the beginning of the period (as in the north-western sub-district).WS trends do not present this problem and makes possible directly analyse the relative variation of the water use within actual rice-cultivated areas.
Figure 10 also shows that in sub-districts A, B and C the decreasing trend is clearly not constant.
The smoothed trends (red dashed lines), show in fact that a marked reduction of WS is evident only WS raw resulted in being quite underestimated with respect to ENR statistics, likely because a single pixel can contain several paddy fields (more than 20, as the area of one pixel is 86 ha, while rice paddies have an average size of 3 ha [52]).Therefore, if paddies included in a single pixel are not flooded simultaneously, FF max may underestimate the fraction of paddies flooded at least once during sowing.In fact, being F j (t) the flooding status (zero for dry, one for flooded) of a single paddy j (j ∈ {1, ..., m}) in a composite image t (t ∈ {1, ..., n}), the real maximum total extent of flooded surface during sowing should be computed as: while FF max was computed as: where ∑ m j=1 F j (t) = FF(t), connecting Equation ( 16) to Equation (3).From this different computation method, it follows that FF max ≤ FF max on all possible conditions.However, FF max is not directly calculable from input 1 × 1 km maps, since single paddies are not discernible; this consideration justifies the use of an empirical correction of WS raw values to balance the observed bias between FF max and FF max .
ENR statistics were essential to obtain unbiased WS measures.Nevertheless, the possibility to estimate this variable from satellite data rather than from farmer declarations is important for several reasons: (i) remote sensing provides impartial and spatialized measures, while a map of the paddies is needed to compute the surface area of declared water-and dry-seeded fields (in this case, the area of cadastral parcels was used, which is often not updated and does not correspond to the actual extent of farming parcels); (ii) remote sensing data allow performing retrospective analysis up to early 2000, hence extending the temporal span of the trend analysis; (iii) estimations from remote data can be automatically updated every year (this will allow further extension of the analysis of the detected changes in the future, without requiring additional field/official data); (iv) the satellite-based estimation can be provided earlier in the season instead of waiting for the collection of farmers declaration, since input MODIS products are freely available two weeks after they are collected by the sensor.In conclusion, satellite data allow providing a near-real-time seasonal view of the water extent and easily broaden the analysis of the flooding condition trends in future years.

Analysis of Changes in the Proportion of Water-Seeded Area
The analysis conducted on the temporal variations of the estimated proportion of Water-Seeded rice (WS) (globally and among sub-districts) provided results quite similar to the ones described for FF avg trends, but with some interesting differences in terms of changes over time.
Results of the ANOVA analysis on the proportion of water-seeded rice fields (Table 6) show that, similarly to what was observed for FF avg , the sub-district is the most important factor in explaining WS variance, but year also has a statistically-significant effect.However, differently from what was observed for FF avg , the interaction between district and year is not statistically significant (F = 1.19, p-value = 0.32).The univariate linear regression over years showed a significant WS decreasing trend of −0.91 ± 0.28% per year (t = −3.27,p-value = 0.0051).This decrease is similar to that of FF avg (−0.86 ± 0.20% per year) and corresponds to a reduction of −14.6 ± 4.5% (−376.8± 1.6 km 2 ) throughout the whole study period.The analysis within sub-districts (see Table 7 and Figure 10) shows that, also in the case of WS, a general decreasing trend is observed in all sub-districts, although some are not statistically significant (i.e., Sub-districts A, B and D).
However, comparing slope coefficients reported in Tables 4 and 7, it can be observed that, while in the case of FF avg , the greater decrease occurred in the northwestern sub-districts, in the case of WS, the situation is the opposite.This apparent contradiction is justified by the fact that the slope of FF avg reduction depends also on the different rice fractional cover.A similar relative reduction of FF avg leads in fact to a larger slope of the regression for sub-districts characterized by a higher FF avg at the beginning of the period (as in the northwestern sub-district).WS trends do not present this problem and make it possible to analyze the relative variation of the water use within actual rice-cultivated areas directly.
Figure 10 also shows that in Sub-districts A, B and C, the decreasing trend is clearly not constant.The smoothed trends (red dashed lines) show in fact that a marked reduction of WS is evident only in the later years (i.e., 2009 onwards), while in the first ten years, no particular trends are evident.On the contrary, Sub-districts E, F and G show constant decreasing trends, as underlined by the almost coincident smoothed and linear trends.These evidences suggest that the lower decrease ratio of WS proportion in the north-western sub-districts is due to a delayed adoption of the dry seeding practice.Further investigations on WS trends during future years are certainly needed to confirm this qualitative finding, and to monitor if the WS proportion reduction will progress in the north-western area with the same intensity already detected for the south-eastern one.
It is also interesting to highlight the anomalous behaviour of 2013 in all sub-districts (more evident in the southern ones).This is due to the atypically strong and continuous spring rain occurred during sowing in that year [79], which forced farmers to sow rice in flooded conditions.This finding is confirmed by ENR statistics (see also Figure 1 of Ranghetti et al. [52]).
Monitoring the trends of FF avg in the next years will be crucial both to i) detect cases of extreme reduction of flooded surfaces in the south-eastern sub-districts (east of the Ticino river), where the flooding extent is already low and a further reduction could lead to a critical shrinking of feeding habitats for herons and egrets; and ii) prevent that the spread of dry seeding in north-western sub-districts reaches the intensities already observed on the south-eastern ones.12)); blue lines represent the linear trends (Table 7); while red dashed lines are smoothed tendencies.Moreover, WS values in the first years of the analyzed period are close to one in Sub-districts A, B and C, while they are already close to 0.5 in the other ones (D, E and F), or even lower (G).
This evidence suggests that the lower decrease ratio of WS proportion in the northwestern sub-districts is due to a delayed adoption of the dry seeding practice.Further investigations on WS trends during future years are certainly needed to confirm this qualitative finding and to monitor if the WS proportion reduction will progress in the northwestern area with the same intensity already detected for the southeastern one.
It is also interesting to highlight the anomalous behavior of 2013 in all sub-districts (more evident in the southern ones).This is due to the atypically strong and continuous spring rain occurred during sowing in that year [79], which forced farmers to sow rice in flooded conditions.This finding is confirmed by ENR statistics (see also Figure 1 of Ranghetti et al. [52]).
Monitoring the trends of FF avg in the next few years will be crucial both to: (i) detect cases of extreme reduction of flooded surfaces in the southeastern sub-districts (east of the Ticino river), where the flooding extent is already low and a further reduction could lead to a critical shrinking of feeding habitats for herons and egrets; and (ii) prevent that the spread of dry seeding in northwestern sub-districts reaches the intensities already observed on the southeastern ones.
Comparing Figure 10 (panel "all") with Figure 1 suggests in fact that the temporal trend of flooded surfaces could be associated with the observed reduction in the population size of waterbirds.The number of nests within rice fields (orange line in Figure 1) showed a decrease between 2000 and 2016 and, similarly to what was observed for WS, which decreased speeded-up since 2010.Nevertheless, in order to understand if there is indeed a causality between the reduction of water-seeded rice area and the number of nests, it would be necessary to study the population dynamics of herons and egrets in combination with other environmental and biological factors.This topic is currently under investigation, and the results of the present study are among the bulk of data under analysis.

Conclusions
This study applied a method to estimate the Flooding Fraction (FF) from MODIS data [52] with the aim to investigate the spatial and temporal trend of flooding within rice fields in the Italian rice district from 2000-2016.Results highlighted that the seasonal averaged flood intensity (FF avg ) and the proportion of water-seeded rice-cultivated area (WS) decreased by −0.86 ± 0.20% and −0.91 ± 0.28% per year, respectively.Moreover, the seasonality of flooding events within the rice-sowing period changed, with a delay in the date of maximum flooding of 0.75 ± 0.31 days per year, suggests a modification in agro-practices and/or cultivated varieties.These changes were found to be widespread throughout the study area, but with strong differences among sub-districts for both FF avg and WS.Analysis of WS trends showed that the spread of dry seeding was already in progress in 2000 in the southeastern parts of the rice district, while it began only recently (about 2009) in the northwestern ones.
During the study period, the rice-cropped area was not reduced, so that the diminution of the standing water areas is to be attributed to the adoption of dry seeding.According to the described results, the adoption of this technique produced a loss of 44% in the surface of paddies where herons and egrets could forage compared to 2000.The loss is greater in the southeastern part of the rice district, reaching up to 76% and 80% in the two eastern sub-districts.Such a reduction may explain the decreasing trend of heron and egret populations observed in the last few years.
These results, and those of concomitant studies on the foraging ecology of the breeding herons and egrets [54,80], will be used in order to assess the effect of environmental and biological factors involved in the recent decrease of these waterbirds.
To limit the massive spread of rice dry seeding, the 2014-2020 Common Agricultural Policy (CAP) already adopted redistributive payments for farmers who implement ecological mitigation (Action 214).This can be done both in Piedmont (Sub-districts A, B, C and D), where Actions 8.2.9.3.2 of the local Programma di Sviluppo Rurale (rural development program) [81] provide support for an anticipated stropping of droughts and for maintaining a partial flooded surface within drought paddies, and in Lombardy (Sub-districts E, F and G), where the attention is focused on the creation of a ditch along each paddy [82] (see Action 10.1.c)to maintain a minimum water reserve without losing efficiency in water management.

Figure 1 .
Figure 1.Total number of nests of the two most abundant breeding waterbirds, the grey heron and the little egret, in north-western Italy during the time window analysed in this study, compared to their predominant foraging habitat[22,23].

Figure 2 .
Figure 2. Localization of the northwestern Italian rice district.The rice-cultivated surface is shown in green, while non-arable land is shown in grey; thick lines and letters denote the sub-districts identified for the analyses (seeTable 1 for details); thin lines denote municipality boundaries.

Figure 4 .
Figure 4. Temporal distribution of DOY(FF max ) over the whole rice district, with smoothed line tendency.

Figure 4 .
Figure 4. Temporal distribution of DOY(FF max ) over the whole rice district, with smoothed line tendency.

272 3 . 2 .
Analysis of the spatial distribution of flooding conditions 273 The spatial distribution of overall flooding conditions during the sowing phase (FF avg ) was 274 not uniform, with higher values detected in the north-western parts of the study area than in the 275 south-eastern one.

277identified.Figure 4 .
Figure 4. Temporal distribution of DOY(FF max ) over the whole rice district, with smoothed line tendency.

Figure 7 .
Figure 7. Distribution of seasonal FF aggregated values (values with NDVI < 0.4) during the entire 2000-2016 study period within each sub-district and the whole district. 280 303

Figure 6 .
Figure 6.Distribution of seasonal FF aggregated values (values with NDVI < 0.4) during the entire 2000-2016 study period within each sub-district and the whole district.

Figure 7 .
Figure 7. Top panels: yearly maps of average flooding presence FF avg maps; bottom panels: 2000-2016 average (left) and coefficient of variation (right) maps.

Version February 16 ,Figure 8 .
Figure 8. Temporal distribution of FF avg pixel values along year and sub-districts (white boxplots) or the whole district (yellow boxplots).Solid blue lines represent the linear regression of equation 11, while red dashed lines are smoothed tendencies, both performed on sub-district based averaged values.

Figure 8 .
Figure 8. Temporal distribution of FF avg pixel values along the year and sub-districts (panels "A" to "G") or the whole district (panel "all").Solid blue lines represent the linear regression of Equation (11), while red dashed lines are smoothed tendencies, both performed on sub-district-based averaged values.

Figure 9 .
Figure 9. Logistic regression between WS ENR and WS raw : calibration points, regression curve with 95% confidence interval, validation metrics and equation.

Figure 9 .
Figure 9. Logistic regression between WS ENR and WS raw : calibration points, regression curve with 95% confidence interval, validation metrics and equation.

Figure 10 .
Figure 10.Temporal distribution of WS values along year and sub-districts.Points are predicted values (equation 12), blue lines represent the linear trends (Table7), while red dashed lines are smoothed tendencies.

Figure 10 .
Figure 10.Temporal distribution of WS values along year and sub-districts.Points are predicted values (Equation (12)); blue lines represent the linear trends (Table7); while red dashed lines are smoothed tendencies.

Table 1
for details); thin lines denote municipality boundaries.

Table 1 .
Characteristics of the rice district and sub-districts."Rice area" is the surface classified as "rice" in the Corine Land Cover images (see Section 2.2.1).
Seasonal FF and NDVI distribution of pixels of the MODIS images (2000-2016) grouped by DOY composite period.NDVI distributions include all the pixels (red dashed line is the proportion of values with NDVI < 0.4), while FF distributions include only selected pixels (NDVI < 0.4).Yellow boxes (grey in the printed version) represent the composites selected for the analysis (see paragraph 2.4.1).

Table 2 .
ANOVA table (type II tests) of the multivariate linear regression (Equation (10)) of DOY(FF max ) over years and sub-districts (with the interaction term).Seasonal FF and NDVI distribution of pixels of the MODIS images (2000-2016) grouped by DOY composite period.NDVI distributions include all the pixels (red dashed line is the proportion of values with NDVI < 0.4), while FF distributions include only selected pixels (NDVI < 0.4).Yellow boxes (grey in the printed version) represent the composites selected for the analysis (see paragraph 2.4.1).

Table 2 .
ANOVA table (type II tests) of the multivariate linear regression (equation 10) of DOY(FF max ) over years and sub-districts (with the interaction term).

Table 1
[52]i) the greater 282 and earlier adoption of dry seeding in the eastern part, leading to lower FF max values.Considering 283 the altitudinal gradient of this area, which decreases from north-west to south-east, it is also evident 284 how the averaged Flooding Fraction decreases from the higher, north-western sub-districts towards 285 the lower, south-eastern ones, in accordance with our knowledge about geographical distribution of 286 water availability in the study area.north-westernsub-districts,beingnearer to rivers, have higher 287 water availability with respect to the south-eastern ones.This probably led to a greater spread of dry 288 seeding in the south-eastern sub-districts, as already pointed out in Ranghetti et al.[52].The map 289 of the coefficient of variation of seasonal FF avg values (Figure6, last panel) also highlights that the 290 interannual variation of standing water presence is higher in the eastern sub-districts, where water is 291 less available.292Figure6 finally highlights that sub-district E (Lomellina) exhibits differences between its 293 north-western (higher flooding fractions), north-eastern (lower values) and southern (intermediate 294 values) parts.This is due to differences in crop fractional cover (rice is dominant in the north-western

Table 3 .
ANOVA table (type II tests) of the multivariate linear regression (Equation (11)) of FF avg over years and sub-districts (with their multiplicative effect).

Table 4 .
Summary of the univariate linear regressions of FF avg over years (each line refers to an independent regression).Intercepts refer to year = 2000 instead of year = 0 to be more intelligible.Stars beside p-values identify significant values (p-value < 0.05).

Table 6 .
ANOVA table (type II tests) of the multivariate linear regression (equation 12) of WS over years and sub-districts (with their multiplicative effect).

Table 6 .
ANOVA table (type II tests) of the multivariate linear regression (Equation (12)) of WS over years and sub-districts (with their multiplicative effect).

Table 7 .
Summary of the univariate linear regressions of WS over years (each line refers to an independent regression).Intercepts refer to year = 2000 instead of year = 0 to be more intelligible.Stars beside p-values identify significant values (p-value < 0.05).