Space-time susceptibility modeling of hydro-morphological processes at the Chinese national scale

1 Hydro-morphological processes (HMP; any process in the spectrum between debris flows 2 and flash floods) threaten human infrastructures and lives; and their effects are only expected 3 to worsen in the context of climate change. One of the ways to limit the potential damage 4 of HMP is to take preventive or remedial actions probabilistically knowing where and how 5 frequently they may occur. The expected information on where and how frequently a given 6 earth surface process may manifest itself is referred to as susceptibility. And, for the whole 7 Chinese territory, a susceptibility model for HMP is currently not available. 8 To address this issue, we propose a yearly space-time model consisting of a Generalized 9 Linear Model of the binomial family. The target variable of such model is the annual pres10 ence/absence information of HMP per catchment across China, from 1985 to 2015. This 11 information has been accessed via the Chinese catalogue of HMP, a data repository the 12 Chinese government has activated in 195X and which is still currently in use. This binary 13 spatio-temporal information is regressed against a set of time-invariant (catchment shape 14 indices and terrain attributes) and time-variant (urban coverage, rainfall, vegetation density 15 and land use) covariates. Furthermore, we include a regression constant for each of the 16 31 years under consideration and also a three-years aggregated information on previously 17 occurred (and not-occurred) HMP. We consider two versions of our modeling approach, an 18 explanatory benchmark where we fit the whole space-time HMP data, including a multiple 19 intercept per year. Furthermore, we also extend this explanatory model into a predictive 20 State Key Laboratory of Resources and Environmental Information Systems, Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences, Beijing, 100101, China University of Chinese Academy of Sciences, Beijing, 100049, China University of Twente, Faculty of Geo-Information Science and Earth Observation (ITC), PO Box 217, Enschede, AE 7500, Netherlands Jiangsu Center for Collaborative Innovation in Geographic Information Resource Development and Application, Nanjing, 210023, China Collaborative Innovation Center of South China Sea Studies, Nanjing, 210093, China Earth Observation Center, German Remote Sensing Data Center, Land Surface Dynamics, Oberpfaffenhofen, 82234, Wessling, Germany Research Center on Flood and Drought Disaster Reduction of the MWR, Beijing, 100038, China State Key Laboratory of Simulation and Regulation of Water Cycle in River Basin, China Institute of Water Resources and Hydropower Research, Beijing 100038, China School of Civil Engineering and Architecture, Southwest Petroleum University, Chengdu, 610500, China one, by considering four temporal cross-validation schemes (Forward-All, Forward-Sequence, 21 Backward-All, and Backward-Sequence), removing the yearly multiple intercept. In the first 22 of 31 temporal replicates, Forward-All is calibrated for 1985 and then used to predict from 23 1986 to 2015. In the second step, a model is calibrated for 1985 and 1986 combined and 24 used to validate the rest of the space-time series. This is replicated up to the last model 25 where the combined data from 1985 to 2014 is calibrated to predict the last year of the 26 HMP presence/absence data. Forward-Sequence also moves in the same temporal direction 27 but the sampling scheme sequentially extracts two years at a time, one for calibration and 28 one for validation. For instance, the first step is trained for 1985 and used to predict 1986; 29 then the second step is trained for 1986 and used to predict 1987. As for Backward-All, and 30 Backward-Sequence, their structure is the same but the temporal direction goes from 2015 31 to 1985. 32 Our explanatory model suggests that the overall number of HMP events per year has 33 increased in the last decade and that the annual susceptibility has subsequently followed 34 the same trend. As for the four cross-validation routines, Forward-Sequence shows excellent 35 performance with an average AUC of 0.83, slightly better than Forward-All, Backward36 All, and Backward-Sequence. From an interpretative standpoint, this implies that the best 37 spatio-temporal prediction we obtained is associated with short-term variations of the HMP 38 distribution and that such variations should be considered in a forward temporal direction. 39 Furthermore, we portrayed the annual susceptibility models into 30 maps, where the 40 south-east of China is shown to exhibit the largest variation in the spatio-temporal proba41 bility of HMP occurrence. Also, we compressed the whole spatio-temporal prediction into 42 three summary maps. These report the mean, maximum and 95% confidence interval of the 43 spatio-temporal susceptibility distribution per catchment, per year. 44 The information we present has a dual value. On the one hand, we provide a platform 45 to interpret environmental effects on HMP at a very large scale, both spatially (the whole 46 Chinese country) and temporally (31 years of records). On the other hand, we provide 47 information on which catchments are more prone to experience a HMP-driven disaster. 48 Hence, a step further would be to select the more susceptible catchment for detailed analysis 49 where physically-based models could be tested to estimate the potentially impacted areas. 50

Hydro-morphological processes (HMP; any process in the spectrum between debris flows and flash floods) threaten human lives and infrastructure; and their effects are only expected to worsen under the influence of climate change.Limiting the potential damage of HMPs by taking preventive or remedial actions requires the probabilistic expectation of where and how frequently these processes may occur.The information on where and how frequently a given earth surface process may manifest can be expressed via susceptibility modeling.For the whole Chinese territory, a susceptibility model for HMP is currently not available.To address this issue, we propose a yearly space-time model built on the basis of a binomial Generalized Linear Model.The target variable of such model is the annual presences/absences of HMP per catchment across China, from 1985 to 2015.This information has been accessed via the Chinese catalogue of HMP, a data repository the Chinese Government has activated in 1950 and which is still currently in use.This binary spatio-temporal information is regressed against a set of time-invariant (catchment shape indices and geomorphic attributes) and time-variant (urban coverage, rainfall, vegetation density and land use) covariates.Furthermore, we include a regression constant for each of the 31 years under consideration and also a three-years aggregated information on previously occurred (and notoccurred) HMP.We consider two versions of our modeling approach, an explanatory benchmark where we fit the whole space-time HMP data, including a multiple intercept per year.Furthermore, we also extend this explanatory model into a predictive one, by considering four temporal cross-validation schemes.As a result, we portrayed the annual susceptibility models into 30 maps, where the south-east of China is shown to exhibit the largest variation in the spatio-temporal probability of HMP occurrence.Also, we compressed the whole spatiotemporal prediction into three summary maps.These report the mean, maximum and 95% confidence interval of the spatio-temporal susceptibility distribution per catchment, per year.The information we present has a dual value.On the one hand, we provide a platform to interpret environmental effects controlling the occurrence of HMP over a very large spatial (the whole Chinese country) and temporal (31 years of records) domain.On the other hand, we provide information on which catchments are more prone to experience a HMP-driven hazard.Hence, a step further would be to select the most susceptible catchments for detailed analysis where physicallybased models could be tested to estimate the potentially impacted areas.For transparency, the results generated in this work are shared in the supplementary material as GIS (geopackage) files.

Introduction
In this work, the term hydro-morphological process (HMP; Wang et al., 2021aWang et al., , 2021b))using was used to address a class of earth surface phenomena where solid and fluid phases of a gravitationally-driven moving mass are not well determined.Thus, this class refers to a broad spectrum of processes spanning from debris flows to debris floods and floods.The reasons behind such initial disclaimer are due to the nature of the dataset we used and further explanations will be provided later in the text.
This class of HMPs includes some of the most frequent and damaging natural hazards, and their occurrence shows a close relationship with climatic changes (Blöschl et al., 2020).HMPs have increasingly been reported to threaten human lives and infrastructure in recent years (e.g., Guo et al., 2018).To prevent or limit losses, it is crucial to estimate where and when these processes may occur.In turn, this enables administrations to plan ahead and mitigate future risks (Rossi et al., 2019).
HMPs are extremely rapid phenomena.Just few hours or even less are needed between the triggering heavy rain and their manifestation (Marchi et al., 2010).They can also be generated by snow-melt but it is generally the intensity and duration of precipitation that control the process through the water discharged over a given area.Then, the overland flows follow the river network, entrains all sorts of debris and leaves it strewn especially when the runoff intersects urban areas (Norbiato et al., 2008).In such cases, roads may be blocked, drainage systems clogged, cars trapped, lives lost and property destroyed (Mahmood et al., 2017.For this reason, HMP prediction models are primarily implemented in a physically-based framework where one can reliably introduce the rainfall input and simulate the process by accounting for topography and soil hydrological characteristics (Tramblay et al., 2010).This is usually performed specifically for small areas (Rozalis et al., 2010) but recent advancement have led to develop similar applications on much wider regions, simulating different types of HMPs from catchment (Javelle et al., 2010) to country-wise scales (Gourley et al., 2017), and even up to continental scales (Paprotny et al., 2017).These different levels of details all share a common structure where a design storm is used as the input.The design storm can be either inferred from long time-series of rainfall data via extreme value statistics (Umer et al., 2022).Or, it can be directly plugged in by using near-real-time rainfall data obtained from meteorological forecasts (Collier, 2007).As for the remaining information, terrain characteristics are commonly derived from global DEM data or from site-specific LiDAR surveys.Besides, soil parameters are required to describe the hydrological characteristics and the associated ability to retain water or to convert it into runoff (Norbiato et al., 2008).This can be obtained via in-situ tests whenever the area is relatively small (Cenci et al., 2016) and from global estimates such as ISRIC, for large scale assessments (Ragettli et al., 2017).These methods have the inherited ability to produce HMP runout estimates, such as total impacted spatial extent, flow heights, kinetic energy, volumes and more, which are crucial information for engineering design and master plans (Li et al., 2019).However, the applicability of physically-based models inevitably suffers from considerable limitations whenever the study target involves continental to global scales ( Van den Bout et al., 2021), with very few exceptions to this rule (Liao et al., 2012).In fact, for large areas, the required input information is typically quite smooth, assuming it is even accessible.And, collecting suitable geotechnical data is difficult if not impossible (Gaume et al., 2009) over large regions.As a result, a complementary branch in the natural hazard community has developed statistically-based models during the last decades.This methods do not offer the same breath of results produced from the physically-based counterpart (e.g., they do not spatially predict runout-impacted areas nor flow-heights, etc.).However, they provide useful information on areas potentially subjected to HMPs, learning from past events from which spatio-temporal projections are made (Gourley et al., 2013).
The present work fits in the second category.Specifically, the Chinese government has recently completed a long lasting project where all the available information on historical HMPs has been collated for the whole Chinese territory.We use the term HMP specifically because the Chinese catalogue reports a wide spectrum of earth surface processes without explicitly attributing a class.This catalogue starts from reports gathered even from ancient China and it covers the period until 2015.
Because of this wide temporal coverage, the data differs in quality across space and time and the Chinese government has decided to use a more general classification, consistent through time.More specifically, the data collated until 1949 is relatively poor and the situation improves substantially from 1950 onward as the current Chinese government was established.Nevertheless, even from 1950 up to 1980, the data may still have some positional issues because the digital system did not exist (Li et al., 2018).The Chinese HMP report system became standardized after the 1980ies, with more available technologies being used to record the location (latitude and longitude), date and time as well as the losses, expressed either in the number of victims or economical value (Guo et al., 2018).In light of these considerations, we subset the Chinese HMP catalogue extracting all the available information from 1985 to 2015.
We note here that since 1985 we also have access to meteorological digital data collected and aggregated daily from the Chinese rain gauge network.
We use this data to build a space-time HMP susceptibility model.A susceptibility model essentially estimates the probability of occurrence of a given natural process within specific mapping and temporal units.Mapping units constitute the spatial structure under which a given study area is subdivided.They can consist of a regular lattice (usually gridcells or rarely hexagons) or they can represent geographic features such as catchments or administrative units (Carrara et al., 1995;Lombardo et al., 2019;Titti et al., 2022).Irrespective of the specific geometry, a mapping unit represents the object upon which a statistical model estimates the probability of occurrence of the target hazard.As for temporal units, they represent the time span upon which the selected model makes a prediction.For physically-based models, this is typically expressed in hours or days whereas for statistical models this may involve a much larger time span.In this work we opted for a catchment partition, having accessed the most updated watershed delineation of China (Shen et al., 2017).As for the temporal partition, we selected an annual unit of time.
As for the method, we chose a binomial Generalized Linear Model (GLM) assuming that the spatio-temporal population of HMPs across China behaves according to a Bernoulli probability distribution.This procedure is quite common and actually represents the most common practice in the geomorphological literature (e.g., Budimir et al., 2015;Tanyas et al., 2021).
The susceptibility to any surface process is not stationary or timeinvariant (Lombardo et al., 2020).It actually varies through time as the environmental conditions change (Jones et al., 2021).For instance, landscape evolution processes may modify the terrain, hence changing the hydrology of a given area.Similarly, road, settlement and urbanization growth can influence surface hydrology (e.g., Tanyas et al., 2022).This is particularly valid in China where the resident population and urbanization rate increased from 10.64% in 1949 to 59.58% in 2018, according to the National Bureau of Statistics of China.This may have changed the distribution of permeable surfaces in favor of concreted and sealed land covers (see, Gong et al., 2019).Also, climate changes may contribute to vary the HMP triggering conditions through space and time, especially because rainfall regimes have become less diluted during wet seasons and they have become more concentrated in narrow time windows.All these contributing/triggering factors can be accounted for in statistical models.For instance, if climate change and accelerated urbanization control the HMP occurrence distribution, then a space-time statistical model should be able to capture their influences and show a potential increase in HMP occurrences in recent years.
We stress here an important topic before closing this Section.The definition of susceptibility corresponds to a purely spatial probability N. Wang et al. term which does not feature the size of the HMP nor its temporal recurrence (Fell et al., 2008).Conversely, the definition of hazard features all these terms at once (Guzzetti et al., 1999).The models we present in this work do not fall in either of the two definitions mentioned above.In fact, we modeled the occurrence of HMPs throughout the Chinese territory but we also modeled their temporal recurrence.As a result, our model output returns a probabilistic estimate which features both spatial and temporal characteristics.Therefore, we do not exactly produce a susceptibility model because our output contains more information than just a traditional susceptibility.On the other hand, our model does not account for the size of HMPs hence, it does not fully satisfy the definition of hazard.What we did, essentially falls in between these two definitions and could be considered a dynamic realization of the susceptibility through time.To simplify our description to the readers, we will simply refer to our output as susceptibility throughout the rest of the manuscript, knowing though that it is not fully correct with respect to the common literature.
The present manuscript is organized as follows: Section 2 introduces the study area and Section 3 describes the material and methodology framework used in susceptibility modeling.This is followed by a detailed description of the model performance and the resulting susceptibility maps in Section 4. Finally, Section 5 discusses the supporting and opposing arguments on this work.And Section 6 summarizes our final remarks.

Study area
China approximately covers the area between latitudes 18 • and 54 • N, and longitudes 73 • and 135 • E. It is characterized by a vast territory and a complex landscape.Based on geomorphological characteristics, China can be divided into six homogeneous regions (Wang et al., 2020): eastern plains (EP), southeastern hills (SEM), southwestern mountains (SWM), north-central plains (NCP), northwestern basins (NWB), and Tibetan Plateau (TP).About two-thirds of China is covered by mountainous areas (Liu et al., 2018a).The southern China consists of hilly and mountainous terrains, while the western and northern China is dominated by plains and basins.The annual rainfall records are strongly controlled by the distance to the coastline and precipitation amounts gradually decrease from the southeast to northwest of China.The eastern plains and southern coasts are severely influenced by the East Asian Summer Monsoon, where most of China's agricultural land and settlements are located.In this context, only the northwest China has a predominantly arid climate and a lower population density.

Hydro-morphological processes in China
As introduced in the previous studies (e.g., Liu et al., 2018b), the Chinese catalogue of HMPs is a digital collection of events, describing a spectrum of phenomena where a fast moving massconsisting of a illdefined proportion of solid and fluidpropagates across the landscape, potentially causing destruction in its path.As a result, the above mentioned spectrum encompasses processes from debris flows (where the solid and liquid phases are almost equally represented) to flash floods (where the fluid phase is much larger than the solid one).HMP records were collected from multiple sources, including Bulletins of Flood and Drought Disasters in China, historical chronicles, together with official documents issued by local governmental departments.Each HMP record in the database contains information on geographic coordinates, date and time.Here, the geographic coordinates indicate the location where the HMPs generated damages.Overall, the Chinese database reports 24,956 HMPs in the time span of 31 years  with a substantially varying concentration across space and time, with the exception of the western arid to semi-arid sector where essentially no events have been recorded (Fig. 1).

Mapping unit
The nature of the Chinese HMP catalogue implies that the various processes may act on different spatial scales.For instance, debris flows usually have a more limited spatial extent, thus slope-to catchmentbased models are the most suitable to represent the physical expression of these phenomena.Conversely, flash floods can travel much longer distances, therefore covering larger geographic scales.In this case, the most suitable models are expressed at scales that range from slope to N. Wang et al. regional ones.Because of this, choosing the most appropriate mapping unit becomes a crucial step to handle the spatio-temporal dimension of the HMP data.We recall here that a mapping unit, in its most basic form, represents the geographic object upon which the landscape is partitioned.In case of relatively small study areas, examples exist where HMPs are modeled along specific streamlined and neighboring areas by adopting a fine squared lattice.This type of resolution and characterization of the HMPs cannot be used in our case, where the size of the Chinese territory would result in billions of grid-cells or data points.Therefore, in case of such large geographic context, a reasonable spatial partition choice could be represented by administrative boundaries, upon which estimating the probability of HMP occurrences.However, the resulting susceptibility model would neglect the hydrology behind the natural process.In fact, administrative boundaries do not necessarily follow streams or catchment divides, where HMP occurrences can be considered independent or nearly-independent from each other.Therefore, a good solution to represent the spatial scale of HMPs, while respecting the hydrological realization of the natural phenomena, is to consider a catchment partition of the Chinese territory.To support the analyses in this work, we selected the most detailed catchment delineation (the 12th level) from the Hydrological data and maps based on Shuttle Elevation Derivatives at multiple Scales (HydroSHEDS database, https://hydrosheds.org),which partitions our study area into 73,587 catchments.The corresponding distribution of catchment sizes is bimodal (see second panel of Fig. 2) and it spans from 0.1 km 2 to 667 km 2 , with average area of 130 km 2 and a 95% confidence intervalmeasured as the difference between the 97.5 and the 2.5 percentiles of the distributionof 231 km 2 .
As any other mapping unit partition used the context of susceptibility modeling, a pre-processing step is required.The presence/absence information of HMPs is to be assigned to each catchment.To do so, we assign a presence (1) and absence (0) label to catchments where at least one HMP record is contained within a specific temporal unit (see details below).

Temporal unit
As much as the mapping unit choice aggregates HMP occurrences over space, whenever a dataset has a temporal connotation one should also choose a temporal unit.A temporal unit is the time interval through which we aggregate HMP occurrences and assign a suitable presence/ absence conditions.In our case, the HMP dataset has a very fine resolution, with date and time available.However, the environmental properties or covariates we will use in the model (see Section 3.4) do not share the same temporal resolution.For instance, rainfall and temperature are available with a daily resolution across China, vegetation cover and urban development are available on a yearly basis while terrain properties do not exhibit any temporal changes.Therefore, choosing a timescale that allows for meaningful interpretation and suitable data is also crucial.In this context, the coarse temporal resolution of the covariates inhibits our ability to build a finely resolved space-time model.And, in any case, choosing a fine temporal resolution would inevitably increase the computational burden.Thus, choosing a reasonable trade-off is required.Due to the characteristics of some crucial covariates, we chose a yearly temporal unit.Such temporal unit implies that we assign a presence (1) and absence (0) label to catchments where at least one HMP record is contained within a year time window.

Covariate set
HMPs are the result of several interplaying factors.These primarily feature: 1. terrain attributes, for they control the path of the overland flows as well as the availability of material to be mobilized and transported; 2. catchment morphology, for it controls the time of concentration and other hydro-dynamic parameters; 3. soil hydrology, for it controls the interaction of the water with the earth surface; 4. precipitation, for it represents the main trigger; 5. temperature, for it controls evapotranspiration and hence the soil moisture; 6. vegetation density, for it can absorb part of the rainfall discharge and interact with soil through the root system; 7. urbanization, for it may change the natural hydrology both because of impermeable surface placed over permeable ones, and because buildings can also reduce the hydraulic section through which HMPs may flow into.
In the context of space-time modeling, these properties need to be considered both in terms of their spatial distribution and temporal evolution.In fact, some properties will be more stationary over time, whereas some will have a much more rapid rate of change.For instance, at the scale of the Chinese territory, soil hydrology can be considered quite stationary within the 31 years under consideration.Conversely, rainfall, vegetation and urbanization might have a much faster spatiotemporal variation.Therefore, certain properties can be introduced as a single realization (or map) whereas other properties should be accounted for their successive temporal realizations (or maps).
The modeling protocol we implemented makes use of both types of covariates, featuring properties that can be safely considered timeinvariant within three decades: terrain and catchment characteristics as well as soil type and climatic zones.And, also by featuring properties that are explicitly time-variant within the same period: climate, vegetation and human activity, as well as antecedent HMPs (this variable will be explained later in the text).
Due to the size of the study area and the temporal connotation of the database, the number of covariate is inevitably large especially because a crucial step consists of aggregating the covariate values in space (at the catchment scale) and time (at the yearly scale).Due to the numerous data sources, the spatial resolution of the covariate set we chose ranges from 90 m (SRTM, https://earthexplorer.usgs.gov) to 8 km (NDVI, http s://climatedataguide.ucar.edu).To summarize the spatial signal of each covariate, we calculated its mean and standard deviation per catchment.In case of stationary covariates, such as terrain attributes, the spatial mean and standard deviation is a sufficient approximation where the mean reflects the main bulk of the pixel distribution per catchment and the standard deviation highlight the associated variability.These values are kept constant through time.As for catchment morphological indices, one single value is computed per catchment and even in this case, the indices are kept constant through time (they are repeated for each of the 31 years).
For covariates that are nonstationary over time (such as rainfall, temperature and vegetation), we compute the spatial mean per catchment as well as the temporal mean and standard deviation in a year (as per Loche et al., 2022).As for the anthropic signal, the percent of urbanized area with respect to the total catchment size is directly calculated on a yearly basis, hence it does not need any spatio-temporal aggregation.To this purpose, we employed the World Settlement Footprint (WSF) Evolution which outlines at 30 m spatial resolution the global settlement growth from 1985 to 2015 on a yearly basis (Marconcini et al., 2020a).The WSF evolution has been generated by exploiting the recently released WSF2015 layer, which maps worldwide the settlement extent for the year 2015 (Marconcini et al., 2020b).In particular, for each pixel denoted as settlement in the WSF2015, a temporal analysis has been performed by means of historical Landsat-5 and Landsat-7 optical satellite imagery to identify when the construction took place.Here, an iterative procedure has been implemented wherestarting backwards from 2015 -training samples for the settlement and non-settlement class are extracted out of the map obtained at time t and Random Forest binary classification has been employed to outline the settlement extent at time t-1.Ultimately, zonal statistics have been computed to determine yearly for each catchment partition the corresponding total amount of settlement area.
As previously mentioned, we also considered antecedent HMPs, calculated over a time window of three years and binarized into presence/absence conditions per catchment.To do so, we tried to capture some residual dependence over space via antecedent HMPs per catchment, then to carry the spatial signal of the HMPs.In fact, within a relatively short time window, we expect the susceptibility to HMPs to be quite spatio-temporally consistent or stationary.In other words, areas that have experienced HMPs in the recent past are more likely to suffer from HMPs in the near future (Samia et al., 2017).Hence, introducing the information of previously occurred HMPs should better inform the model of this short-term spatial dependence and improve its overall prediction capacity (Lombardo et al., 2020).
A simpler overview of the covariates we considered is provided in Appendix 7. To these reference set of variables, we also add a multipleintercept per year exclusively for a first set of analyses where we will model the whole spatio-temporal data at once (explanation provided below).The inclusion of a yearly intercept is meant to capture residual temporal dependence in the data, under the assumption that climate change may have led to a larger number of HMP occurrences in recent years.

Susceptibility modeling
In this work, because of the vast study area and the long time series, we opted to create a susceptibility model that can feature spatiotemporal characteristics.We do so by considering two types of models, an explanatory one and a set of predictive ones.The explanatory model is a model built by using the whole available information.In our case, it is a model where the entirety of China is taken into consideration together with its 31 years observations.In such a way, one can build a model that can be used for interpretation, to understand the statistical role of every environmental factor with respect to HMP occurrences.However, such models do not have a predictive connotation because no new or unknown data is used to test the classification performance.
We stress here that the natural hazard communityat least the part of it using statistical modelsusually performs calibration by randomly subsetting a percentage of the data over space and test the validation performance over the complementary cases.However, prediction or forecast are terms usually referred to estimates of future occurrences, hence in time.Rarely, studies dedicated to susceptibility models are validated in time (or chrono-validated) (Lombardo and Tanyas, 2020;Cama et al., 2015), mostly because of the inherited complexity of obtaining accurate multi-temporal inventories (Guzzetti et al., 2012).
Because our dataset spans over such a large time window, we actually have the chance to test whether it is possible to forecast future occurrences.Thus, in addition to a traditional 10-fold cross-validation run over the whole data, we have opted to assess the predictive capacity of future HMP occurrences by considering four temporal cross-validation schemes: N. Wang et al.Notably, each of these validation schemes inevitably produces 30 testing outputs, whereas the explanatory model only produces one training output.

Generalized linear models
The vast majority of statistically-based susceptibility models are carried out by using Generalized Linear Models (Budimir et al., 2015;Reichenbach et al., 2018).This class of models assumes that the response variable follows an exponential family distribution such as Gaussian, Poisson, Bernoulli and more.Among those, the Bernoulli case, also referred to as Binary Logistic Regression, corresponds to a model where the target variable can take on only two values.Therefore, a binomial GLM estimates the probability that a given mapping unit belongs to one of the two classes (by standard, this is the class 1, or the class conveying the presence of HMPs, rather than 0).More specifically, a binomial GLM can be denoted as follows: where, the target variable Y is assumed to be Binomial with a probability π of a given catchment to experience a HMP.The β 0 term is the global intercept and β n are the regression coefficients estimated for X n covariates.The logit, or the natural logarithm of the odds, allows for the conversion of the odds into probabilities.This framework allows for continuous and discrete covariates.Each class of a discrete covariate is modeled independently from the other classes, or technically it is assumed to be independent and identically distributed (iid).More specifically, the model will assign a different regression constant to each class separately from the others.Notably, in this work we make use of iid covariates for a multiple yearly intercept for the explanatory reference model.The remaining covariates are all continuous in nature and used as linear properties both in the explanatory and predictive models.

Estimates of confidence intervals
In statistics, any model should allow for inference on a distribution of estimates rather than a single estimate.In other words, obtaining a mean prediction is as important as measuring the uncertainty around that mean value.Therefore, in this work we sought to retrieve both the mean behavior of every regression coefficient and performance metric as well as their estimated variability.
To do so, we present two schemes, one for the explanatory model and one for the validation routines (MOD1 to MOD4).When implementing the explanatory model (we recall here that it is fitted using the whole available information), we have also added a bootstrap simulation step (Efron and Tibshirani, 1994).This step essentially re-samples with replacement the whole dataset and re-fits the same model structure to the simulated dataset.We do this over 100 bootstrap replicates to estimate the sampling distribution of each parameter we store during the explanatory analyses.Besides, we implement the 10-fold cross validation to evaluate the overall performance on the whole dataset.As for the validation routines in MOD1 to MOD4, the variability of the tests is summarized via the 30 estimates, one for each of the 30 years under consideration.

Model evaluation
The primary tool to assess the performance of our HMP susceptibility model consists of the Receiver Operating Characteristic curves (ROC, Hosmer and Lemeshow, 2000) and their integral or Area Under the Curve (AUC, Hosmer and Lemeshow, 2000).The former is the most common threshold independent metric used in classification problems (Rahmati et al., 2019).It is constructed by slicing the probability spectrum at various cutoff, and by computing the confusion matrix at each step.As a result, it is possible to calculate the False Positive Rate or FPR (FP / [FP + TN]) and the True Positive Rate or TPR (TP / [TP + FN]) for each cutoff.The integral of the curve defined by the FPR and TPR pairs calculates from different cutoffs can be then used as an index of performance.Specifically, AUC = 1 indicate a perfect classification, 0.9 < AUC < 1 refers to outstanding performance, 0.8 < AUC < 0.9 marks excellent performance whereas 0.7 < AUC < 0.8 are acceptable results.Any AUC value from 0.7 to 0.5 indicates a range of poor performance down to results comparable to a random classification.
We make use of the AUC throughout the manuscript.We also implement a Jackknife test in the validation scheme (Lombardo et al., 2016a).A Jackknife test is essentially divided into two steps.The first one runs single (j th ) variable models whereas the second runs all-butone-variable (j − 1) models.In both cases, the AUC is calculated to offer a comprehensive summary of covariates contributions.Single variable models highlight stand-alone performance of specific covariates in explaining HMP occurrences.All-but-one-variable models highlight performance drop resulting from the removal of one single covariate at a time, with respect to a full model using them all at once.Notably, the validation scheme in this work includes training and testing 30 temporal models.As a result, we have run 30 Jackknife tests, one for each year from 1985 to 2015.

Explanatory model and its cross-validation
In this section, we reported the regression coefficients obtained from a susceptibility model built by using all the available spatio-temporal information.These estimates were used to interpret the relation between HMP occurrences and environmental conditions (or covariates).Firstly, each regression coefficient is characterized by a distribution of values which have been retrieved from 100 nonparametric bootstrap replicates.Fig. 3 summarizes each model component.Among the continuous covariates (see Fig. 3), climatic indices (e.g.RAIN_Tσ_Sμ, AnnualRAIN_Sμ), terrain attributes (e.g.PLC_σ, SLP_σ), catchment morphology (e.g.form factor) present notable positive regression coefficients.In addition, catchments located in Central temperate and South temperate zones also suffer more from the HMPs.More details on the interpretation of this covariate effects will be provided in Section 5.

N. Wang et al.
Besides, we made use of an iid effect for each year, whose result is shown in Fig. 4. The year-specific regression constants show an interesting pattern.For each year from 2002 to 2014, all regression coefficients are significantly positive and the whole distribution is quite distant from the zero line (between 0.5 and 1) with an exception of 2004.As for each year in the period between 1985 and 2001, the regression constants are also estimated with a positive median coefficient, although some of them appear to be not significant (the distribution of regression constants also show negative values).Besides, the regression coefficients vary around zero.More details on the interpretation of this temporal iid effect will be provided in Section 5.
To complete the analyses on the whole spatio-temporal domain, we also run a 10-fold cross validation.We recall here that a 10-fold cross validation implies randomly partitioning the whole data population into ten complementary subsets, each time extracting 90% and 10% for calibration and validation, respectively.Fig. 5 presents the performance of the 10-fold cross-validation scheme.Specifically, panel 5a reports 10 ROC curves obtained by using 90% of the spatiotemporal HMP data; and panel 5b reports ROC curves obtained by testing over 10% of the spatiotemporal HMP data.The respective mean AUC values do not significantly change, as they both returned 0.84.This attest both for excellent goodness-of-fit and prediction-skill according to Hosmer and  2000) as well as a indicating robust results with differences that can be distinguished only at the third decimal place.

Temporal validation routines
Here we present the four temporal validation schemes described in Section 3.5.For each temporal validation scheme, we summarize the model performance in Fig. 6.All models are reported with a mean temporal AUC greater than 0.82.We recall that this value corresponds to excellent performance according to the AUC classification system proposed by Hosmer and Lemeshow (2000).However, two distinct patterns arise in the four temporal validation routines.The AUCs obtained for each year in MOD1 and MOD3 appear quite smooth.In MOD1, this is also associated with a downward shift in AUC when comparing calibration and validation performances (Fig. 6a).As for MOD3, calibration and validation performance largely overlap, with the exception of the period in between 2009 and 2015 where the validation routine shows a significant drop in predictive capacity (Fig. 6c).In case of MOD2 and MOD4, the AUC values estimated for each year present a much rougher temporal variation.Between these two validation schemes, MOD4 less accurately predicts the HMPs in the last years of our AUC time series (Fig. 6d).As for MOD2, a similar difference in performance between calibration and validation is shown for the initial years of our HMP time series (Fig. 6b).However, the initial years from 1986 to 1989 contain less HMP occurrences, thus a relatively low performance in this period is much more acceptable than a relatively low performance in the latest years.In light of these considerations, and because of a slightly better performance overall, we consider MOD2 (or Forward-Sequence) as the best validation scheme compared to the other three.
We stress again that a close look at MOD2 in Fig. 6b highlights some fluctuations in the AUC time series for the validation whereas the calibration appears much more stable through time in terms of estimated performance.This is better presented in Fig. 7 where we show 30 ROC curves, one for each year.The panel (a) corresponds to the training ROC curves and aside for a few years, they consistently overlap through time.As for the validation shown in panel (b) a marked spread can be seen in the curves spanning from 1986 to 2015.We note that the relatively poorer performance registered at the start and end of the time series also     N. Wang et al. correspond to two years with a relatively lower number of observations.Conversely, the other relatively low AUC values between the two endpoints always appear in the following year containing very large numbers of HMP occurrences.This may be due to the fact that an abrupt increase in HMPs, may induce some variations in the estimated dependence between HMPs and covariates.In turn, this may also induce variations in the susceptibility patterns, which may end up not matching the HMPs of the subsequent year (likely to be representative or closer to the average Chinese susceptibility pattern).This explanation fits well the year 1998.That year was characterized by an exceptionally large number of HMPs in China, and the temporal validation of 1999 returned the poorest performance we observed across the whole temporal sequence.Notably, such temporal variations in performance has been similarly shown in other studies, where the authors reported that effect of climate change may be responsible for large uncertainties in the prediction of HMPs (e.g., Collier, 2007).
To provide a comprehensive overview of the model structure and covariates' role in the temporal validation, we performed a suite of Jackknife tests (Jiao et al., 2019).We recall here that Jackknife tests are essentially replicates of a reference model whose structure is perturbed by either building single-variable (only-one-variable) models for each of the covariate in the reference structure.Or, by removing one covariate at a time (all-but-one-variable) from the whole set of covariates.Many example of Jackknife tests exist in the susceptibility literature, but they have been limited so far to a pure spatial domain (see for instance, Lombardo et al., 2016b;Ramos-Bernal et al., 2019).Here, because we consider both spatial and temporal dimensions, we iterated the only-onevariable and all-but-one-variable models thirty times, once per year from 1985 to 2015.
Fig. 8a presents the AUC obtained via only-one-variable models.It indicates that most of the terrain attributes, climatic indices, and antecedent HMPs could contribute to a model with an AUC greater than 0.6.At the same time, the all-but-one-variable models in Fig. 8b indicates that removing either of SLP_σ, form factor, elongation ratio, RAIN_Tσ_Sμ, and antecedent HMPs from the model will induce an obvious AUC drop.

Regression coefficients
In addition to assessing model performance, another crucial step in any modeling protocol is to evaluate how reasonable regression coefficients may be from an interpretative standpoint.In this work, we already summarized a similar information for our benchmark fit.Nevertheless, regression coefficients estimated for the temporal validation scheme could shed light on the variability that each covariate effect may exhibit through time.Here, we assigned the yellow colour for a positive β value, which indicates the probability of HMP occurrence will increase by a factor equal to the exponential of the β value.Conversely, the blue colour indicates a decrease.
Among the terrain attributes, the standard deviation of slope (SLP_σ) and plan curvature (PLC_σ) play an important role in controlling the estimated probability of HMP occurrences (Fig. 9).In terms of catchment morphology, form factor and elongation ratio show a positive effect.Most soil types present non significant and negligible contributions to the thirty cross validation schemes, with the exception of Sandy Clay which appears to be negatively correlated to HMPs, although with a slight negative influence.Furthermore, catchments located in Central temperate, South temperate, and Central subtropics zones appear to be more prone to HMPs than the others.
The summary presented above reports the role of time-invariant properties.As for time-variant covariates, AnnualRain_Sμ showed the largest significant and positive effect out of all the climatic indices, followed by RAIN_Tσ_Sμ (the intra-annual rainfall variance within a given catchment).The 3-years antecedent HMPs in a given catchment also appeared to be significant and to increase the susceptibility estimates.
Notably, the summary of the covariates' effects shown above is quite static as it overlooks the temporal variation that each model component exhibit as the temporal-validation is performed.To complement this information, in Fig. 10 we show the temporal evolution of the regression coefficients belonging to covariates that appeared to be significant in Fig. 9.
More specifically, to better distinguish the variance of the covariates' Fig. 9. Regression coefficients obtained by MOD2.
N. Wang et al. effects through time, we split Fig. 10 in two panels, according to the magnitude of the regression coefficients.Panel (a) summarizes β coefficients whose magnitude through time ranges from − 0.5 to 0.5, whereas panel (b) presents the same information for β coefficients whose magnitude through time ranges from − 2.5 to 5. Most of the covariates in both panels indicated a constantly similar effect on HMP occurrence, whereas, few covariates showed a large variation through time.For instance, the annual rainfall (AnnualRAIN_Sμ) indicated an increasing positive influence from 1985 to 2014.However, the variance of NDVI (NDVI_Tσ_Sμ) within each year showed a decreasing effect before 1990, after which the trend flattened until the end of the time series.We stress here that accounting for the time-varying signal of the vegetation, here through the NDVI, is rarely performed in the landslide literature, with the exception of few cases (e.g., Schmaltz et al., 2017) Nevertheless, the covariates which exhibited the largest influence and variation through time all correspond to climatic indices, especially those associated with rainfall (see AnnualRAIN_Sμ and RAIN_Tσ_Sμ in Fig. 10b).

Susceptibility mapping
HMPs susceptibility maps generated via MOD2 are drawn in Fig. 11 from 1996 to 2015.These have been classified into very low (VL), low (L), low to medium (LM), medium to high (MH), high (H), and very high (VH) according to break points that have been set as the 2.5%, 25%, 50%, 75%, and 97.5% percentiles of the whole probability spectrum combined.In other words, to reclassify the numerical susceptibility into classes, we have concatenated all the space-time HMP probability estimates into a single vector, from which five percentiles have been extracted to ensure a common colour palette among the 30 maps.
Looking at the time series of susceptibility maps (Fig. 11), at the beginning of the study period probabilities are generally lower, especially in the western sector.Besides, as the time series evolves towards recent years, the probability spectrum appears to shift towards higher susceptibility classes.More specifically, catchments with very low probabilities of HMP occurrences mainly appear from 1986 to 1988; whereas catchments presenting very high probability of HMP occurrences characterize the south-east sector of China since 1997.
To summarize the space-time susceptibility information depicted in  Looking at the susceptibility patterns depicted in the mean and maximum maps, two macro-regions stand out the most.The western sector of China appears to be consistently estimated as non susceptible.There, the susceptibility appears to be generally confined within the first 10% of the national distribution.On the contrary, the south-eastern sector appears to be generally the most susceptible.There, most of the catchments are associated with susceptibilities estimated above 70% of the national probability distribution.Notably, few exceptions exist to this observation due to the existence of large plains, where catchments are generally gentler in topography.Other catchments with high HMP susceptibility, albeit lower than the south-east, can still be found in central, north-east and southern China.
Interestingly, the 95% confidence interval mapwe recall here to be computed as the difference between the 97.5th and 2.5th percentiles of the spatio-temporal probability spectrum shown in Fig. 11 marks analogous geographic patterns to the mean and maximum maps.This is an indication of the robustness of our model.In fact, this means that areas with low susceptibilities are estimated with similar values through time.Conversely, areas with high susceptibility exhibit a much larger degree of variation through time, as expected just by looking at the raw data in Fig. 1.The last panel of Fig. 12 depicts seven cluster drawn from the maximum susceptibility in the same fig..These have been manually interpreted on the basis of expert opinion.Clusters I to V correspond to regions are affected by monsoon.The reason to split I and II are due to the difference of terrain and annual rainfall whereas the reason to split I and III into two parts is due to the mountain range that acts as a strong topographic barrier.More specifically: • Cluster I: the region with the most severe erosion due to the topographic control; • Cluster II: the region mostly affected by monsoon; • Cluster III: less annual rainfall, Loess Plateau affected by widespread gully incisions; • Cluster IV: this sector of China shows a relatively large proneness towards HMP although the rainfall intensity due to the incoming monsoons in this area is much lower than the precipitation discharged to the south.This is primarily due to the local rough topography which contributes to increase the probability of HMP occurrence; • Cluster V: plains with widespread flat terrains; • Cluster VI: distinct characteristics attributable to the Taklamakan Desert and the Tarim Basin; • Cluster VII: sparsely populated area corresponding to the Changtang Plateau and Qinghai Hoh Xil Plateau.
Fig. 12 is meant to compress the spatio-temporal susceptibility information in the geographic space.To do the same for the temporal case, we went back to Fig. 11 and computed for each year the areas assigned with one of the six susceptibility classes.From these, we generated a stacked barplot (see Fig. 13) reporting the proportion of China associated with one of the six classes, showing the evolution through time from 1986 to 2015.What stands out the most is that the areal percentage of catchments with very low (VL) susceptibility decreased sharply in the first three years.This effect is mostly due to the fact that as the time series progressed, more HMP have been recorded, which generally leads to a higher probability of HMP.On the opposite side of the probability spectrum, the proportion of China associated with very high (VH) HMP susceptibility can be seen to have increased from 1997 onward.We remind here the reader that despite these changes may appear small in a simple graphical summary such as Fig. 13, in reality a variation of even just 1% of the total Chinese territory involves several hundreds thousands of km 2 and several hundreds actual catchments.

Supporting arguments
This work estimates and investigates the spatio-temporal variation of HMP susceptibility patterns over China.Because of the vast space-time domain, many options exists on how to build and validate a spacetime predictive model (Lombardo et al., 2020).
We chose a binomial GLM, which we calibrated and validated through different strategies.The first strategy we used exploited the whole space-time domain, from which catchments with high variations in slope steepness and planar curvature appear to increase the overall susceptibility to HMPs.The influence of slope with respect to HMPs is widely acknowledged in the literature.However, for analyses expressed at the catchment scale, the effect of the terrain steepness may be lost.This may be the reason why in our model, the positive role of the slope steepness is expressed in terms of standard deviation, a common proxy for topographic roughness.A similar reasoning can be inferred for the standard deviation of the planar curvature.
Unsurprisingly, another positive contribution is carried by the rainfall patterns, expressed through the RAIN_Tσ_Sμ and the AnnualRAIN_Sμ (see Fig. 3).It should be noted that the spatio-temporal rainfall signal is carried in the model via four summary statistics of the precipitation over the mapping (catchment) and over the temporal (year) units.This is certainly the reason behind the overall negative contribution estimated for RAIN_Tμ_Sμ.In fact, in any multivariate analysis, whenever slightly dependent covariates interact with each other, the estimation of their sign and amplitude can also depend on each other presence within the model.Because the four rainfall aggregation measures stem from the same spatio-temporal information, it is safe to assume that some degree of dependence can exist among the four we computed.Therefore, the overall influence of rainfall on HMP occurrences should be interpreted as the combined effect of the four covariates and their estimated regression coefficients, which returns an overall increasing effect of the HMP susceptibility as the rainfall increases (see Fig. 3 and note the following median values: β RAIN_Tμ_Sμ = − 0.75, β RAIN_Tσ_Sμ = 0.80, β RAIN_TA_SA = − 0.04, β AnnualRAIN_Sμ = 0.67).
As for the temperature, the effect is much clearer there, as all the three summary statistics derived from the original spatio-temporal temperature signal appear to have a negative contribution to the model.This means that at increasing temperatures the probability of HMP occurrences consistently decreases in space and time, irrespective of the three components at hand.
We also stress here the relevance of antecedent 3-years HMPs.This idea stems from the fact that certain types of hazard persist or cluster both in time and space, hence by featuring antecedent occurrences in the model can help predicting future HMPs.This concept is not new in landslide studies (Samia et al., 2017), although a similar approach has not been tested yet when modeling HMPs statistically.
An additional and equally interesting contribution to the model was carried by human interference.Other researchers have already pointed out a similar consideration (Plate, 2002), which we tested in this work by including the presence of build-up area per catchment and per year (Marconcini et al., 2020b).The expansion of human settlements has a dual effect in our model.On the one hand, it undeniably constitutes an interference with the environment, potentially leading to HMP occurrences (Duncan, 2013).On the other hand, human expansion also means that a larger number of people are being exposed to hazards (Cutter et al., 2018).This in turns may bring some degree of bias in the HMP inventory because events that occur in non-inhabited areas may not be recorded, especially due to the size of the study area.Conversely, events that occur in inhabited areas are much more likely to be recorded.
As regards the temporal validation scheme we tested, it is important to justify why we chose MOD2 as our best and further presented it to the readers.When looking at performance, not only the central tendency (mean or median) but also the variation around it constitute a relevant criterion.The variation is essentially described as the difference between the highest and lowest performance.Among the two terms, we chose the lowest performance, together with the median AUC, to be our primary mean of selecting the best temporal validation scheme.In fact, in an ideal situation one should avoid selecting a model that can poorly perform even as rare as this may occur.Therefore, MOD2 has become our best validation scheme for it both provides a median value comparable to MOD1, MOD3 and MOD4.And, it returned a minimum AUC much higher than the other temporal validation routines.
In terms of covariates' influence on HMP susceptibility, MOD2 offers a slightly different perspective than the first exploratory tests.The morphological characteristics of the catchments largely contribute to the HMP susceptibility (see form factor and elongation ratio in Fig. 10.And even more interestingly, RAIN_T σ _S μ and AnnualRAIN_S μ not only dominate the probability estimates to a much larger extent than any other covariate.But, they also show a quite distinctive increasing trend through time.
An increasing trend through time is also shown for the probability of HMP occurrence at the scale of the whole Chinese territory.This becomes evident when looking at Fig. 11.This illustration showcases early years where the western Chinese sector was characterized by very low susceptibilities, whereas the same territories have shifted towards low susceptibilities in few decades.Similarly, the south-eastern landscape N. Wang et al. initially appeared to be dominated by low to medium susceptibility with fee high hotspots.However, in the thirty years under consideration, these areas have shifted to predominantly highly susceptible, with localized very high hotspots.Fig. 13, summarizes the same information, offering an even more straightforward graphical intuition.The very low susceptibility areas have rapidly decayed in surface extent and from the early 2000's, have become barely visible.In a opposite fashion, the high and very high susceptible areas have rapidly increased and nowadays characterize a relevant number of catchments.
Ultimately, we decided to remind the reader that the susceptibility we estimated for the whole Chinese territory is actually much finer in resolution than what it looks like in the previous fig.s(e.g., Zhao et al., 2018) (The smallest catchment area can be 0.1km 2 in this study).To do this, we have selected eight important and large catchments.In Fig. 14 we show the HMP susceptibility estimated for each of those catchments via MOD2, and specifically through the maximum probability of HMP occurrence shown in Fig. 12.Looking at Fig. 14 becomes evident that our model is built on a very detailed spatial partition of the Chinese landscape.And, within each of the eight selected major catchments, it is possible to further distinguish susceptible sub-catchments that upon which local administrations can made decisions to ensure the safety of local inhabitants.

Opposing arguments
The model we present is both spatial and temporal in nature.Among the suitable space-time models we have chosen a relatively simple one, a binomial GLM.Instead of this, we could have opted for a binomial Generalized Additive Model (GAM, Titti et al., 2021) extension to account for possible nonlinear covariates' effects.And, to include potential variables acting at a latent level, hence requiring complex CAR (conditional autoregressive) or SPDE (stochastic partial differential equation) components to be featured as well.We maintain that our choice has proven to be a valid option, for both our spatio-temporal cross-validation and temporal validation schemes returned AUC values well above 0.8, the threshold for excellent binary classifiers according to Hosmer and Lemeshow (2000).
However, even in this case, one could say that the performance shown in Fig. 6 appear to be quite constant irrespective of the temporal validation scheme we opted for.In this sense, the catchment and year units we selected may have played a role.In fact, using catchments to model HMP occurrences may not always be the ideal situation.These processes may act at a scale smaller than the catchments' extent.Thus, like any other relatively coarse mapping unit, a catchment may act as a spatial smoother.The same consideration could be made for the N. Wang et al. temporal unit.The triggers of HMP may act at a daily scale or at best at a monthly scale, if we consider the effect of antecedent rainfall.But, a yearly model, could act as a temporal smoother.Thus, a higher resolution both in the spatial and temporal dimensions could lead to a larger and more interesting variability in model performance and related interpretations.This being said, the model we propose is built with the most reasonable mapping unit if we consider the continental extent of the Chinese territory.And, it is built with the most reasonable temporal unit if we consider the 30-years epxression of the catalogue.
Some may also argue that on a 30-year long record, the accuracy of the inventory may have drastically changed in recent times.As a result, the inventory may be biased towards a larger number of HMP records at the end of the time series.We maintain that the HMP inventory is reliable and should not be affected by this type of bias, at least to the best of our knowledge.In fact, the Chinese government has supported the initiative of creating this inventory long before the 1980'ies.And, by the starting year of our time series (1985), the recording protocol had already been standardized at the whole Chinese territory scale.Surely, we cannot entirely disregard the possibility of some sort of bias due to the size of the study area.We know for a fact (and already shared the information with the readers in Section 4) that the western sector of the Chinese territory is devoid of large settlements.This may imply that the lack of HMP record in the region and the subsequent low susceptibility estimated there (see Fig.s 11 and 12) could be due (to some extent) to a lack of interest rather than a real absence of HMPs.We have tried to investigate this potential issue, by checking local news and other source of information.But, we have not found records of HMPs in the region.Thus, we can only assume the inventory to be reliable.
Furthermore, researchers have studied the effects of incomplete inventories onto the modeling results.And, they confirmed that ignoring flaws in the mapping procedure can lead to geomorphologically incoherent results.Nevertheless, the propagation of these flaws onto the model estimates can be limited or even adjusted by adopting specific counter-bias modeling strategies (e.g., Steger et al., 2016).These strategies have been already extensively described in the literature and for additional details we suggest reading Steger et al. (2021).Here, for a brief summary, we stress that counter-bias strategies boil down to including the bias into the model when it comes to the fitting phases but they require to remove potential covariates, sensitive to the bias, when it comes to the prediction phase.This is referred to as a zeroing-out procedure of the regression coefficients estimated for specific categorical properties, again, sensitive to the bias (see Lin et al., 2021, for a specific application in China).Alternatively, and this is also an option we are currently evaluating, one could opt for deep learning approaches.In recent years, these have been reported to produce stable results irrespective of different sampling strategies and potential biases in the data (Wu et al., 2017) and especially in high-dimensional data such as ours (Li and Wang, 2020).
Nevertheless, this should not be our case, or at least, we do not have indications that the widespread distribution of catchments with HMPs' absences is due to a human bias or if it is due to natural causes.In fact, the spatio-temporal dataset we have built features a much larger number of catchments without HMP records rather than catchments recorded with actual HMPs.In other words, we used an unbalanced presence/ absence dataset.In turn, this affects the estimated probability spectrum, resulting in a positively skewed susceptibility distribution.We maintain that this is a realistic pattern for such a vast spatio-temporal domain where probabilities of HMP are generally very low, with the exception of few catchments that are very susceptible (see Frattini et al., 2010).However, it is worth mentioning that the geomorphological community often operates with a balanced dataset of presence and absence cases.This in turn makes the probability spectrum much more normallydistributed and centered at around 0.5 (Rossi et al., 2010).Both approaches are valid, although creating a balanced presence/absence dataset distorts the global intercept (Lombardo and Mai, 2018) making the interpretation of the probabilities valid only in a relative sense rather than the common notion of probability available in any other statistical application (Petschko et al., 2014).Therefore, we maintain that our unbalanced choice is valid and suitable to treat such a large spatiotemporal domain.

Conclusions
The Chinese territory hosts a vast and diverse landscape that in the last thirty years has been struck by thousands of hydro-morphological processes.Such processes, spanning from debris flows to debris floods and floods have been monitored and recorded in a multi-temporal digital archive thanks to a Chinese program centrally coordinated but enacted by local administrations.In this work, we explore and exploit this archive to produce the first catchment-scale-based HMP susceptibility model of China, from 1985 to 2015.
We distinguished seven macro-regions where the estimated probability of HMP occurrence can be interpreted and explained.The southeastern regions are the most susceptible to HMPs, primarily because of the monsoon control on precipitation regimes.This observation of a spatial patterns is strictly entwined with the temporal observation that the susceptibility estimates tend to increase in recent years.This may be due to the fact that climatic changes have narrowed the duration of storms and increased their intensity.
For transparency, we are sharing the file containing the susceptibility estimates for each year under consideration.Although we cannot share the raw data, the method we propose is reproducible in any other geographic context.For this reason, we consider our work an example of continental-decadal scale HMP susceptibility.We stress here that other examples currently present in the literature have all been built by using a grid-cell partition of space, where each grid-cell has a resolution in the order of kilometers.Therefore, their actual use is hindered by the fact that over several kilometers, the landscape can contextually feature floodplains as well as mountain ridges.
The strength of the model we propose is due to the temporal breath of data we examined on a yearly basis.This characteristic, makes it unique together with the representation of the geographic space at the catchment scale.As a result, the information we provide is expressed at a scale that respects the geomorphology and hydrology of the phenomena under consideration.
The operational advantage of the model we propose is that it be considered as a first order indicator of catchments under threat.Geotechnical or more generally physically-based models aimed at simulating HMP genesis, propagation and deposition are extremely expensive because of their data requirements.However, one can use our output to select catchments with large probability of HMP occurrence, where to prioritize data acquisition for engineering solutions related to HMP occurrences.susceptibility modeling of hydro-morphological processes at the Chinese national scale", by Nan Wang, Weiming Cheng, Mattia Marconcini, Felix Bachofer, Changjun Liu, Junnan Xiong and Luigi Lombardo has no conflict of interest of any kind.

Fig. 2 .
Fig. 2. Probability density distribution of catchments sizes in China, computed from the most detailed catchment delineation (the 12th level) from HydroSHEDS database (https://hydrosheds.org, last access on 7th December 2020).Red dots correspond to HMP occurrence locations.(For interpretation of the references to colour in this fig.legend, the reader is referred to the web version of this article.)

Fig. 3 .Fig. 4 .
Fig. 3. Regression coefficients estimated through the explanatory model built by using the whole HMP spatio-temporal information across China.The covariates shown in this fig.are continuous in nature.The red dash line corresponds to zero or no-contribution to the model.Boxplots shown in blue indicate a median negative correlation to HMPs while yellow indicates a median positive one.(For interpretation of the references to colour in this fig.legend, the reader is referred to the web version of this article.)

Fig. 6 .
Fig. 6.Each panel corresponds to one of the four temporal validations we tested.The line plots report the AUC time series from 1985 to 2015.The boxplots summarize the AUC variation over the thirty years.
N.Wang et al.

Fig. 8 .
Fig. 8. Jackknife test for covariates.Only-one variable models are shown in the left panel and all-but-one variable models in the right panel.Red line indicates the corresponding mean value of all combinations.Blue boxplots indicate a covariate-specific median AUC lower than the mean AUC computed for all covariates.Yellow boxplots correspond to higher covariate-specific median AUC.(For interpretation of the references to colour in this fig.legend, the reader is referred to the web version of this article.)

Fig. 11 ,
Fig. 11, we further generated three maps aimed at delivering the mean and the maximum susceptibility together with the 95% confidence interval (see Figs. 12a, 12b and 12c respectively).Looking at the susceptibility patterns depicted in the mean and maximum maps, two macro-regions stand out the most.The western sector of China appears to be consistently estimated as non susceptible.There, the susceptibility appears to be generally confined within the first 10% of the national distribution.On the contrary, the south-eastern sector appears to be generally the most susceptible.There, most of the catchments are associated with susceptibilities estimated above 70% of the national probability distribution.Notably, few exceptions exist to this observation due to the existence of large plains, where catchments
1. Forward-All or MOD1: This validation procedure starts by calibrating our binomial GLM (more details in Section 3.5.1)overaspecificyear(e.g., 1985)and testing over the remaining time series(e.g.,  1986-2015).In the second step, the previous reference year is combined with the next(e.g., 1985 and 1986)to predict HMPs in the remaining years(1987 to 2015).This moving window moves one year at a time until completion of the time series.2.Forward-Sequence or MOD2: This validation scheme iteratively calibrates over a specificyear (e.g., 1985)and predicts only the next (e.g., 1986).In the second step, the calibration aggregates the subsequent year (e.g.,1985 and 1986)and predicts only the next (e.g., 1987).This is repeated until the completion of the time series in 2015.3. Backward-All or MOD3: This validation scheme is analogous to MOD1 but it is implemented in the opposite temporal direction.Specifically, we calibrate over the last year (2015) and predict the whole time series backward (from 1985 to 2014).In the next step we then calibrate aggregating the information of the previous year (e.g., 2015 and 2014) to predict the remaining time series (1985 to 2013).This operation is backwardly repeated until the completion of the time series in 1985.4. Backward-Sequence or MOD4: this model is analogous to MOD2 but again it is implemented in the opposite temporal direction.This means that the calibration starts in 2015 and it is used to predict the previous year only (2014).Then the calibration integrates the information from the previous year (2015 and 2014) to predict only one step back in time (2013).This operation is repeated backwardly until the time series is completed in 1985.