A study of land subsidence problems by ALPRIFT for vulnerability indexing and risk indexing and treating subjectivity by strategy at two levels

A decline in water table is widely known to induce subsidence due to a subsequent loss of hydrostatic forces in the form of land settlements including both gradual downward movements and sudden sinking without significant horizontal movements. Subsidence is a highly variable process and has a wide range of impacts including adverse effects on safety, interruptions of services and transportation, as well as environmental impacts and damage. Research on subsidence is topical and some recent works use experimental techniques to understand the complexity of the problem. For instance, Li et al. (2014) reports on a set of experiments on the deformation behaviour of sands under abstraction–recharge cycles. Ma et al. (2010) identify the hysteresis effect or time-dependent lag in deformation when the water level recovers. Wei et al. (2017) investigate experimentally the deformation of samples taken from a confined aquifer and aquitard with lithology of silty clay, silt and fine sands and they characterise elastic–plastic deformation for all samples. Studies of this nature contribute to the better understanding of subsidence and complement earlier investigations, e.g., Bouwer (1977), who shows through simple calculations that a decline of 10 m in water table can induce as much as 5–50 cm subsidence and hence this may be taken as a rule-of-thumb. The diversity of topical research on subsidence is outlined below but the primary aim of the paper is to develop risk management tools for land subsidence when data availability is sparse.

Existing research trends are summarised in Table 1 in broad terms, which categorise them as: (i) mathematical models use modelling capabilities of groundwater flows to predict the decline in water table and thereby to facilitate studies on subsidence and their impacts (e.g., Modoni et al. 2013); (ii) remote sensing techniques are used to monitor land subsidence to investigate effects of hydrogeological settings (e.g., Gao et al. 2016); and (iii) use GIS techniques to study subsidence problems. The latter group is suitable for cases with sparse data but Nadiri et al. (2018a) argue that there is not much published work yet. They introduce ALPRIFT as a new technique, where past techniques are often simplistic, e.g., the hazard risk assessment index techniques (Huang et al. 2012). The paper precludes the techniques in the first and second categories but is focused on innovating ALPRIFT, as detailed below.

The paper builds on vulnerability mapping by the ALPRIFT framework as presented by Nadiri et al. (2018a) and Sadeghfam et al. (2020a), and tests different strategies of risk mapping problem, where mapping produces relative values suitable to cases with sparse data using largely general-purpose data not collected for such studies. The success of the concept of risk is often indicated by the diversification of its application and adaptation of its definitions in a variety of ways. Generally, some of the approaches to calculate risk include: (i) probabilistic approach which quantifies risk by multiplying the consequences of hazard with their likelihood, as reviewed by Khatibi (2011); (ii) combination of load and resistance, as reviewed by Tung et al. (2005) and Sadeghfam et al. (2018); and (iii) indexing risk as a product of vulnerability and hazard (as reviewed in Living with risk (2004)). The probabilistic techniques providing absolute values of risk are often referred to as risk quantification techniques, which use the probabilistic approaches, but these are precluded in the paper for their extensive data requirements. For sparse data, risk indexing (mapping) techniques are used and the paper transforms vulnerability indexing by ALPRIFT into a risk indexing problem.

ALPRIFT is the acronym of seven data layers, which comprise: aquifer media (A), land use (L), pumping of groundwater, recharge (R), aquifer thickness impact (I), fault distance (F) and decline of water table (T). A further data layer comprises results of Interferometric Synthetic Aperture Radar (InSAR). Sadeghfam et al. (2020a) examine closely these data layers and note that the ALRIF data layers (five of seven) have the capacity only to prevail at local scales in accounting for geogenic and natural processes; whereas the PT data layers (the remaining two) prevail at the system-wide scale owing to the hydraulic properties of water acting as a conduit for the processes stemming from anthropogenic activities. For instance, the data layer A is bound to vary at a local scale without any correlation between its occurrence throughout a study area; whereas the P data layer has the capacity to impact system-wide under the right conditions through hydraulic processes. Notably, the paper refers to ALRIF as passive vulnerability index (PVI), whereas PT is referred to as active vulnerability index (AVI).

The distinction on passive and active processes inherent in ALPRIFT (i.e., PVI and AVI) is a justification for developing a risk indexing capability but this needs to be tested under different conditions. The concepts underlying PVI and AVI lead to formulating two problems: (i) the additive model of PVI (ALRIF data layers) and AVI (PT data layers) to serve as the total vulnerability indexing (TVI) problem; (ii) the multiplicative PVI and AVI model is developed as a novel risk index (RI) problem. A critical examination of these new indices is vital to ensure that PVI, AVI, TVI and RI are properties contained in ALPRIFT data layers but not noise. The paper formulates a modelling strategy at two levels. At Level 1, genetic algorithm (GA) is used to optimise the weights of data layers without testing its performance. At Level 2, supervised artificial neural network (ANN) is used, in which the target values are ascertained through data fusion to serve the training and testing phases. The data fusion approach presented by the paper is also and is based on catastrophe theorem.

Besides building evidence for PVI, AVI, TVI and RI that represent significant signals, further novelties in the paper are as follows. (i) A modelling strategy is formulated at two levels but with different selections as normally carried out by the author's ongoing programme of research, e.g., vulnerability indexing of land subsidence (Nadiri et al. 2018a, 2020; Sadeghfam et al. 2020a), vulnerability to groundwater contamination (Nadiri et al. 2018b), groundwater level prediction (Moazamnia et al. 2019). (ii) The data fusion in the paper applies to combining PVI, AVI or RI with InSAR results to produce more consistent results, the algorithm for which is through the innovative use of catastrophe theory, deriving its basis from the decision theory given by Cheng et al. (1996). Data fusion is used widely and the best example is the human mind capable of combining and making sense from different senses. Notably, artificial intelligence (AI) techniques (GA, ANN, fuzzy logic) are now established modelling tools and therefore the paper suffices just to their specification. The strategy is investigated at Tasuj aquifer, which suffers from decline of water table and consequently land subsidence.

This section presents the rationale for transforming ALPRIFT, presented by Nadiri et al. (2018a), into a risk mapping technique, which comprises a new thinking on ALPRIFT.

ALPRIFT is treated through a new thinking, the rationale for which comprises: (i) a breakdown of the ALPRIFT data layers into PVI (ALRIF data layers) and AVI (PT data layers); (ii) the additive model of PVI and AVI provides the TVI problem and their product the RI problem.

The appropriate general-purpose and site-specific data for each of the seven data layers at each pixel are processed at the preliminary stage towards forming datasets, in which the minimum possible SVI is 23 and its maximum possible value is 230. The best practice procedure for the data is given later.

The domain of the influence of each of the ALPRIFT data layers in Equation (1) can be local or system-wide. Similar to Sadeghfam et al. (2020a), the paper argues that a local domain of influence would refer to the correlation among a number of pixels as a measure of PVI of the ALRIF data layers. Their common properties include: (i) prescribed ratings for ALRIF are such that higher SVI values are indicative of a greater potential for subsidence; (ii) they are related to natural geogenic processes; and (iii) while it is reasonable to expect correlation at local scales, owing to a diversity of interplaying factors, system-wide correlations are unlikely. The remaining two data layers of PT are related to the action of the water table, in which P is an anthropogenic factor and T is a direct response to it. These together have a potential for system-wide domain of effects and account for activating hydrostatic pressure through withdrawing from any control volume. Withdrawing hydrostatic pressure triggers active subsidence processes and hence AVI, but can also be referred to as active load or pressure. Notably, the R data layer has also the potential for system-wide effects but it is included in ALRIF, as it accounts for natural processes.

The likelihood dimension concerns the loss of hydrostatic pressure and amplifies the potential for compaction and consolidation towards subsidence. This is accounted for by P (Pumpage) and T (water table). These effects are acted spontaneously to trigger subsidence.

Although vulnerability or risk index is not dimensionally equivalent to subsidence measured in length, each of these quantities are normalised in the range of 0–1 so that they can be mapped onto one another. In this way, the measured subsidence, either in situ or by satellite, is treated as proxy for either vulnerability or risk. Sentinel-1 satellite SAR data are observed information records on Earth's surface, which go back to 2010. They are processed by InSAR using phase difference between two SAR observations. Phase is defined as the fraction of one complete sine wave cycle and is often highly correlated with terrain topography and therefore SAR images together with the phase information contain amplitude information related to the strength of radar responses.

The interferometric wide (IW) swath products in this study comprise a 250 km swath with a spatial resolution of 5 m × 20 m, which include three sub-swaths captured by terrain observation with progressive scans SAR (TOPSAR). This technique provides homogeneous image quality throughout the swath with uniform signal-to-noise ratio (SNR) and distributed target ambiguity ratio (DTAR). The procedures in InSAR processing lead to producing surface elevations in eight steps and are detailed in Figure 1.

The paper uses GA, ANN and data fusion together with the cusp type of catastrophe functions and therefore they are specified in this sub-section for serving as building blocks for the formulation of a modelling strategy on PVI, AVI, LVI and RI and these are below.

The paper uses GA, originally given by Holland (1992), as an optimisation technique to identify the values of the weights. It generates high-quality solutions to optimisation and search problems and is a widely used bio-inspired optimisation method by emulating such genetic operators as mutation, crossover and selection (Mitchell-Olds 1996). Its search process uses four evolutionary processes: (i) initialisation, (ii) selection, (iii) crossover, (iv) mutation (Davis 1991). The paper uses a best practice procedure, in which the initial population evolves as follows: (i) the population in each iteration forms a generation with an assessment of the fitness of each individual in the population; (ii) the genome of each individual evolves and forms a new generation through a modification (recombined and mutated randomly); (iii) the new generation of candidate solutions is used in the next iteration of the algorithm; and (iv) the process is terminated at the maximum number of generations or when the modelled fitness level reaches a designated value. The formulated objective function to identify the weights is given later.

ANNs emulate the working processes in the brain and are parallel information processing systems through a set of neurons for transforming input signals to output signals by using a system of weights expressed as activation functions. The neurons in a layer are not connected to each other but each one in a layer is connected to all the neurons of the next layer. The paper uses a multi-layer feed-forward perceptron (MLP) network, where its topology comprises an input layer, a hidden layer and an output layer. Learning the weight parameters between the input–hidden and hidden–output layers uses the Levenberg–Marquardt (LM) algorithm, hence MLP-LM. The topology of the hidden layer is identified through a trial-and-error procedure. Implementations of ANN are detailed by Haykin (1999) and ASCE (2000).

The modelling strategy formulated by the paper innovates by introducing a data fusion technique that uses catastrophe theory rooted in the decision theory by Cheng et al. (1996). Data fusion is defined by See & Abrahart (2001) as the amalgamation of information from different data sources, which is widely practised in electrical engineering. It is characterised at signal level, pixel level, feature level and decision level (Jiang et al. 2011), each with its own definition and procedures. Abdelgawad & Bayoumi (2012) present another categorisation as low-level fusion, medium-level fusion, high-level fusion and multilevel fusion.

Data fusion models are formulated for inference, estimation, feature maps, aggregation, abstract sensors, classification and compression (Abdelgawad & Bayoumi 2012). A variety of techniques supports data fusion and they include statistical matching, Bayesian inference (as reviewed by Chen & Han 2016), moving average filters (Rodriguez et al. 2015) and grey relational analysis (Bai et al. 2016). The paper innovates by using cusp catastrophe to combine results of a lower level model for PVI, AVI or RI by pairing them with InSAR results as inputs. The outputs of PVI, AVI or RI serve as targets for the Level 2 model.

While the paper uses the scoring system in ALPRIFT, it treats the inherent subjectivity in weights by introducing a modelling strategy at two levels. They comprise: Level 0: pre-process ALPRIFT data layers and prepare the InSAR layer; Level 1: learn the weights associated with each data layer using GA and prepare PVI, AVI, TVI or RI maps; and Level 2: carry out data fusion by using PVI, AVI or RI and InSAR results to produce target data for ANN models; train and test ANNs using ALPRIFT data layers but under direct supervision.

The paper refers to the activities prior to the implementation of the modelling strategy as Level 0 and these comprise any set of activities prior to the risk/vulnerability mapping tasks. Activities at Level 0 are concerned with addressing data availability; deciding the scope of the study and decisions made on transforming raw data into datasets, as well as formulating a modelling strategy.

The modelling strategy at Level 2 is rather complex and comprises the following steps: (i) data fusion for the production of target quantities of PVI, AVI and RI; (ii) formulating three supervised ANN models; and (iii) using the three ANN models for mapping PVI, AVI and RI. These are outlined below and illustrated in Figure 2.

It is noted that each of Equations (6a)–(6c) has an output, which use the normalised value of PVI, AVI, RI and InSAR at each pixel but are adjusted through the n^th root. The choice of the 2nd root for a and the 3rd root for b signify one scheme but other schemes are worth exploring.

Step 2 – Supervised ANN models: Three ANN models are constructed for mapping PVI, AVI and RI. The inputs to these three ANN models use the rated data layers and the target values are fed from Step 1 at Level 2. The models are then trained and tested through the usual trial-and-error procedure to identify their appropriate topology.

Step 3 – Mapping: The maps of the three ANN models are produced for PVI, AVI and RI but TVI mapping is carried out as the sum of PVI and AVI.

The performances of each mapping activity for PVI, AVI, TVI and RI (ALPRIFT, GA and three ANN models) are studied by the receiver operating characteristic (ROC) curve and the area under curve (AUC). Swets (1988) developed this matric to measure the accuracy of a diagnostic system. It divides the events regarding the diagnosis into four groups as true positive (TP), true negative (TN), false positive (FP) and false negative (FN). In the subsidence problem, an event refers to hazard, vulnerability and risk indices, while diagnosis refers to Earth's displacement. The ROC curve plots FP proportion versus TP proportion at various threshold settings. A preferable ROC curve is the one with deviation towards the upper left corner. AUC quantifies the accuracy of the ROC curve, defined as the proportion of the area beneath the ROC curve to the total area. Thus, AUC varies between 0.5 and 1. AUC = 0.5 indicates no discrimination or equal FP and TP proportion while its ROC curve lies along the major diagonal line. On the other hand, AUC = 1 indicates perfect discrimination and its ROC curve follows the upper and left axes.

The performances of the ANN models for the mapping of PVI, AVI, TVI and RI use determination coefficient (R²) and root mean square error (RMSE) metrics. The value of R² for ‘perfect’ model fits is 1 but when its value approaches zero, its performance become poorer as an indication of no correlation. The best RMSE value is 0 but has no upper limit, although the lower the value the better the performance.

The actual values of PVI, AVI, TVI and RI at each pixel are variables but this hampers extracting the patterns for an understanding of their inherent variabilities. Therefore, the results are communicated by their classifications into bands using the Jenks optimisation method. It is a data clustering method to identify the optimum arrangement for different classes through finding the minimum average deviation from the class mean by maximising the deviation of each class from the means of other classes. It reduces the variance within the classes and maximises the variance between classes as the main goals of the Jenks optimisation method (Jenks 1967). This paper uses five bands as follows: Band 1 (low), Band 2 (relatively low), Band 3 (moderate), Band 4 (relatively high), Band 5 (high). Using this technique reduces the subjectivity associated with indices’ classification by expert judgement.

Tasuj (also known as Tesi) plain is located in the west of East Azerbaijan province, northwest Iran and forms the northern coast of Lake Urmia (see Figure 3(a)), currently distressed due to water level decline over recent years since 2000. The plain is bounded by Mishov and Qatir Uchan Mountains at its north, Shabestar plain at its east and Salmas plain at its west. The aquifer and plain formed from alluvial deposits and cover an area of 295 km² (see Figure 3(b) and 3(c)). In general, the diversity of the formations in the basin is wide and all three groups of rocks including igneous, metamorphosed and sedimentary are observable in the geological formations from Precambrian to the recent era. Figure 3(c) also illustrates the spatial location of fault in the south-east of the study area.

This aquifer is unconfined and serves as the main source for domestic and agricultural uses in the city of Tasuj and the villages within the plain. The aquifer and seasonal watercourses flow to Lake Urmia in the direction of north-to-south (see Figure 3(a) and 3(b)). Notably, the general slope of the plain decreases from the north-east, Qatir Uchan Mountains, to south-west, shoreline of Lake Urmia. There are 69 qanats, 41 springs and 144 deep and semi-deep tube wells which discharge an estimated annual volume of 16 × 10⁶ m³ from the aquifer (East Azerbaijan Water Authority (EAWA) 2012). According to the Emberger (1930) method, the study area is characterised as an arid and cold climate zone, and the maximum and minimum temperature are, respectively, 33 °C and −11 °C at Tasuj station in a ten-year period (EAWA 2008–2018). The average annual precipitation is 290 mm during a ten-year period (provided by Meteorological Organization of Iran 2008–2018).

The decline slopes are 2.3 and 8.6 cm per month for Lake Urmia and Tasuj aquifer, respectively, using linear regression analysis. The concern over the rapid decline of groundwater levels at the aquifer is attributable to policy failure with serious environmental and ecological impacts, e.g., saltwater intrusion and land subsidence.

Best practice procedure for the seven ALPRIFT data layers is summarised in Table 2 along with the specification of raw data sources and required processing, similar to Nadiri et al. (2018a). The preparation of data layers L and T (land use and water table) in the paper differs from Nadiri et al. (2018a), as they used the decline in one-year period, whereas the paper uses the trend of water table decline as a more representative measure of decline.

Land use was prepared by image processing through Envi software in three steps, pre-processing, processing and post-processing. The pre-processing step applies geometric and atmospheric corrections by the following steps: (i) use normalized difference vegetation index (NDVI) to identify different land uses; (ii) define the determinable land uses as representative of local knowledge and image interpretations by the supervised classification approach to classify different land uses based on the maximum likelihood method. The results are evaluated using local knowledge and expert judgement with further modifications, if necessary.

Table 4 represents six InSAR single look complex (SLC) images from Sentinel-1A satellite acquired during April 2017 to April 2018. These images were processed five times as illustrated in the last two columns. The final result of InSAR processing is cumulative displacement for five interferometry analyses as discussed in the section ‘Processing InSAR’.

The square grid sizes selected in this study are 20 m × 20 m and this renders the total number of pixels to be 736,740. The dataset for the three ANN models were divided randomly into the training and testing phases in the proportion of 80 and 20%.

Four sets of results are presented: (i) at Level 0, the results comprise the individual rated values (see Figure 4), as well as the processed InSAR results; (ii) at Level 1, the GA results are presented for PVI, AVI, TVI and RI; (iii) at Level 2, the ANN results are presented for PVI, AVI, TVI and RI; and (iv) these results are compared.

The seven ALPRIFT data layers were prepared by processing the available raw data as summarised in Table 2 based on the flowchart in Figure 1. The rated output data layers are presented in Figure 4 for the study area, as well as Figure 5(a), which displays land subsidence in cm/annum for the period from 2017.06.21 to 2018.04.05. Attention is drawn to Figure 4(f) and 4(g), which reveal the variations in P and T data layer with hydraulically connected strips at their overlaps with expected high AVI, as discussed in due course.

The mapping results are classified into five bands by using the Jenks optimisation method (Jenks 1967), in which Band 5 represents high vulnerability/risk, and Band 1, that of low vulnerability/risk. A visual comparison of ALPRIFT in Figure 5(b) with InSAR results in Figure 5(a) indicates areas of strong convergences but also strong divergences, particularly at the strips of PT displaying some hydraulic connectivity. Based on the performance metrics of ROC/AUC in Figure 6 with AUC = 0.56, the performance of ALPRIFT is fit-for-purpose and good enough to serve as the basis for benchmarking the improvements at Levels 1 and 2, although improvements are needed to reduce inherent noise for defensible decision-making. Level 1 modelling results are presented in the next section.

The results at Level 1 use GA, the default parameters for which are given in Table 5 together with the identified weight values for each data layer. The comparison in the table is not directly possible, as the weights of ALPRIFT are assumed to change in the range of 1–5 but those of GA between 1 and 7; notably, ALPRIFT data layers are quite divergent and this may be a reflection of the subjectivity in the scoring system.

GA is used to reduce inherent noise for a defensible decision-making and the results presented in Figure 5(c1)–5(c4) for values of PV, AVI, TVI and RI at Level 1, respectively, are compared visually with both Figure 5(a) (InSAR) and 5(b) (the ALPRIFT benchmark). Even a visual comparison indicates significant convergences particularly along strips of hydraulic connectivity due to PT data layers and some divergences.

The convergence between Figure 5(c2) (AVI values) and 5(c4) (RI values) and those of Figure 5(a) (InSAR) are evident; as well as those between Figure 5(c2) (AVI values) and 5(c4) (RI values) with Figure 4(f) (pumpage). The visual evidence for the role of water to act in strips of the study area to actuate risk is evident and so is the successful role of data fusion. A visual comparison of the results in Figure 5(c2) (AVI values) and 5(c4) (RI values) shows that AVI strongly influences RI along a certain area but in the remaining parts the concordance between PVI and TVI is stronger. Arguably, existence of AVI is a significant learning emerging from the methodology presented by the paper and the role of data fusion is demonstrable.

The curvature of the ROC curve and the performance metrics of AUC in Figure 6 shows improvements in the AUC value for PVI, AVI, TVI and RI to make them defensible in decision-making. The results provide significant evidence for existence of PVI and AVI.

The results at Level 2 are presented in Table 5 using ANN, for which the topology of the three ANN models were decided through a trial-and-error procedure. It also presents the performance metrics of the selected topology for PVI, AVI and RI, according to which there is a slight deterioration in the quality of the performance metric from the training phase to the testing phase, but this is expected. Overall, the models are fit-for-purpose in terms of RMSE and R² performance metrics for studying the performances of PVI, AVI, TVI and RI.

The ANN mapping results, classified into five bands using the Jenks optimisation method (similar to Level 1 models) are presented in Figure 5(d1)–5(d4) for values of PVI, AVI, TVI and RI at Level 2, respectively, and compared visually with both Figure 5(a) (InSAR) and 5(b) (ALPRIFT), as well as with Figure 5(c1)–5(c4). They show notably more convergences and more reduction in divergences. In particular, the convergences between Figure 5(d2) and 5(d4) and those of Figure 5(c2) and 5(c4) are evident, as well as those between Figure 5(d2) and 5(d4) with Figure 4(f). The visual evidence for the role of water to act in strips of the study area as the agent of actuating risk is evident and this is reinforced by the curvature of ROC and the performance metrics of AUC in Figure 6, which shows improvements in the AUC value for PVI, AVI, TVI and RI to make them defensible in decision-making. A visual comparison of the results in Figure 5(d2) (AVI) and 5(d4) (RI) further shows that AVI at Level 2 strongly influences RI at Level 2 along the strip but in the remaining parts, the concordance between PVI and TVI is stronger, as per Figure 5(d1) and 5(d3). The results provide significant evidence for existence of PVI and AVI.

The overview of Figure 6 is that it captures the emerging narrative on the signals inherent in ALPRIFT and the two levels of modelling to improve on the signals as follows. (i) Level 0: the signals in basic ALPRIFT, despite containing subjectivity in rates and weights, are significant as per AUC value of 0.56 and its concave curvature but better information extraction techniques are needed for more defensible decision-making and hence the modelling strategy at two levels. (ii) Level 1: GA renders AUC values of 0.61, 0.63, 0.67, 0.67 for PVI, AVI, TVI and RI, respectively; and the curvatures of ROC become more concaved. (iii) At Level 2, ANNs render AUC values of 0.68, 0.66, 0.70 and 0.71 for PVI, AVI, TVI and RI, respectively; and the curvatures of ROCs become even more concaved. Hence, the modelling strategy makes a difference in a defensible way to extract inherent correlations within the ALPRIFT data layers.

The proportions of the areas swept by Band 1–Band 5 are delineated and presented for each set of results and displayed in Figure 7. The results evidently show that the worst case of Band 5 is currently significant as per each set of results but increased learning tends to unravel a greater proportion for adverse bands. Therefore, searching for more improved modelling strategy is important. The majority of the study area seems to be located in Bands 1–2, but an ongoing investigation by the authors (yet to be published) shows that the situation is worsening.

The convergences and divergences between SVI using ALPRIFT and InSAR results are also the focus of attention, but the reasons for the divergence are not clear-cut. However, it is noted that the correspondence between PVI and TVI, as well as that between AVI and RI are obvious. Notably, the central strip along the north–south direction acting as a hydraulic corridor to actuate risk is in evidence. A further patch at the north-west seems to be high in AVI. The rest of the areas are rather low in PVI and TVI and therefore low in AVI and RI.

The performance metrics including ROC/AUC signify convergences and divergences between InSAR results and PVI, AVI, TVI and RI values. A further measure is ideal to provide a visual aid to understand the convergences and divergences and these are presented in terms of the scatter diagrams for the ANN results, as well as those for the residuals. The results are presented in Figure 8(a1)–8(a3) for the scatter diagram and Figure 8(b1)–8(b3) for the scatter diagram of the residuals. The three scatter diagrams (Figure 8(a1)–8(a3)) display similar shapes, approximately elongated triangular in the range of 0.2–0.7. Values in the vicinity of the 45° line mark better performances but those bulging out of the line display significant divergences that require some explanations. Likewise, three residual scatter diagrams (Figure 8(b1)–8(b3)) display three similar shapes and in the range of 0.2 to −0.3. Values in the vicinity of the line of 0.0 residuals mark better performances but those bulging out of the vicinity show high divergence normally where InSAR values are high. The figures for TVI do not follow the above patterns.

Although the improvements from Level 0 to Level 2 are notable, as per results in Figures 6 and 8, there is room for improvement and one to achieve this is through testing different strategies, as well as more improved measured data. The pixels at low bands are liable for more divergences. The exception for the TVI shapes are likely to stem from not developing an ANN model but deducing it from PVI and AVI.

The field evidence for land subsidence is the basis for ground truthing and a glimpse of the situation on the ground is displayed in Figure 9. As the vicinities of observation wells are often filled by granular material (gravel) after subsidence, the reported values of subsidence from field observation do not provide sufficiently reliable data to test the developed methodology. Additionally, there are other uncertainties associated with the material of the piezometers.

The programme of research by the authors in the paper is towards developing a GIS-based ‘working tool’ for mapping subsidence at plains over aquifers. Arguably, the current status of the reported research on ALPRIFT, PVI, AVI, TVI and RI related to subsidence is a strong proof-of-concept, as classified by NASA (see: https://www.nasa.gov/sites/default/files/trl.png) corresponding to TRL4 or higher. The future direction of the authors’ work is to raise these techniques to Technical Readiness Level 9 (TRL9).

The rates and weights of seven ALPRIFT data layers are prescribed that were suggested by Nadiri et al. (2018a). These values were set to cover different hydrogeological settings of study areas. However, the values are inevitably subjective as they are set by expert judgement and this is a normal, as acknowledged by Nadiri et al. (2018a) and the present study. This explains why the problem of enhancing model accuracy by learning from observed values is topical for all frameworks, such as DRASTIC and ALPRIFT is no exception. The present study aims to reduce inherent subjectivities by a modelling strategy. Evidently, the ALPRIFT framework is generic and applies to all sites; model enhancement strategies are also generic and can be applied to other sites; but the identified rates and weights are site-specific because they are learned from a set of site-specific data.

In addition to the generic nature of the ALPRIFT framework, it is in its infancy and the study areas investigated are yet to expand. It is necessary to investigate more areas with complex hydrological settings such as confined or multiple aquifers. In these cases, although prescribed weights and rates are still generic, some data layers should be prepared in different ways. For example, water table decline is substituted by the decline in piezometric level, or where aquifer crops out should be considered to prepare the recharge data layer.

The authors are exploring different ways of dividing ALPRIFT data layers into PVI and AVI and testing various modelling strategies. Their ongoing research, both reported in the paper and by Sadeghfam et al. (2020a), confirms that the inherent signals are significant in the general purpose ALPFRIFT data layers for having AUC values in the range of 0.6–0.7 to support existence of PVI, AVI, TVI and RI as new quantities. They are currently investigating ways of improving the strength of the identified signals by testing different modelling strategies in two layers. The basis of the techniques developed in the paper may be captured by two terms, heuristics and framework. The term heuristics refers to any broader techniques where their rational justification is judged by the efficacy of the results to enable modellers to discover or learn about data to gain a new insight but normally without any theoretical or empirical basis. The particular choices made in the paper are all on a heuristic basis and without any theoretical or empirical justifications. Also, the term framework signifies the cases, where there is no theory, empirical knowledge or even measurements to justify the selected data layers.

PVI (ALRIF) extracts information on land subsidence with local significance and AVI (PT) on that of anthropogenic activities with system-wide significance. Although ROC/AUC indicates that signals are significant, other combinations of data layers are also feasible, e.g., ALIF and RPT with respect to the role of water creating a system-wide impact.

The characterisation of the result as fit-for-purpose and defensible is discussed in detail by Khatibi et al. (2020) to invoke criticality that modelling results are not the end but means to an end and therefore the underlying assumptions have to be examined, as follows. At Level 0, the authors had directly observed subsidence values but the observation wells were not associated with GPS facilities and therefore the observation values were found unsuitable for ground truthing. At Level 1, the paper uses GA for the optimisation problem as it is widely used, but other techniques should also be investigated. At Level 2, data fusion was investigated by one possible scheme without sensitivity studies, but other techniques should also be investigated. Khatibi et al. (2020) refer to modelling strategies of this nature as ‘inclusivity in multiple models’ (IMM) and argue that there is a mathematical condition to ensure the improvements in the model accuracy. There are a wide range of strategies possible and this itself is the subject of research. The choice for the modelling strategy presented in the paper was a demonstration of the feasibility without maximising the extracted information from the data.

Overall, with the benefit of hindsight, the treatment of TVI in terms of a simple addition of PVI and AVI seems a poor choice simply based on the scatter diagram and the residuals displayed in Figure 8(a3) (scatter diagram for RI) and 8(b3) (scatter diagram of residuals for RI). This does not undermine the modelling results produced in the paper but just that TVI values at Level 2 can be better. The authors’ future plans include testing the performances of learning from inclusive multiple models with different data fusion techniques.

The paper presents evidence on the proof-of-concept for transforming SVI into a RI capability for plains over aquifers, particularly for where there is sparse data availability including uncertain or imprecision inherent in data. The methodology is formed by collective contribution of ALPRIFT data layers using a scoring system of prescribed rates and prescribed weights. The paper divides the data layers into two groups: (i) the five ALRIF data layers account for PVI measuring local resistance to variations; (ii) the two PT data layers account for AVI measuring system-wide effects; (iii) the additive model of PVI and AVI renders the TVI capability equivalent to basic ALPRIFT; and (iv) the multiplicative model of PVI and AVI renders the RI capability.

The methodology is implemented through an IMM strategy at two levels, as follows: (i) at Level 0, preliminary decisions are made; (ii) at Level 1, GA is used as an optimisation technique to identify the weight values for vulnerability and risk mapping of the study area; and (iii) at Level 2, a three-step modelling procedure is devised, which comprises: Step 1: the data fusion of vulnerability and risk indices with InSAR results produce outputs that are fairly placed to serve as target values; Step 2: supervised ANN models are constructed for vulnerability and risk mapping, which use the rated ALPRIFT data layers together with the weights produced at Level 1 and target values as per Step 1 and thereby to produce a trained and tested models; Step 3: the models are applied to produce maps of the study area for vulnerability and risk indices.

A study of the results provides further evidence that in spite of subjectivity in ALPRIFT, it is capable of capturing significant signals from site-specific data but they need to be improved for a defensible decision-making. The two-level modelling strategy formulated by the paper provides evidence that the splitting of the ALPRIFT data layers into two groups for mapping PVI and AVI produce improved signals in vulnerability and risk mapping. As such, there are strips of the study area where the effects of PT are transmitted beyond local pixels and where risk mapping is more reliable; whereas there are areas, where active PT effects are not significant and remain limited to local impacts.

Data cannot be made publicly available; readers should contact the corresponding author for details.