Periodico bimestrale
Anno XIX, numero 88
Sett./Ottobre
ISSN 1128-3874
Ambiente e territorio

Validazione e performance di carte di suscettibilità all’instabilità dei versanti ottenute con metodi deterministici. Il ruolo dei modelli di infiltrazione. Il caso di Antrodoco (Rieti).

Valerio Vitale – Istituto Nazionale di Statistica Vittorio Chiessi - Istituto Superiore per la Protezione e la Ricerca Ambientale

La valutazione e la mappatura della propensione al dissesto costituisce un obiettivo necessario al miglioramento della sicurezza di persone e dei vari elementi antropici che compongono la realtà territoriale. L’espansione di attività economiche e di insediamenti residenziali che prevedono significativi investimenti necessitano di strumenti di valutazione che minimizzino gli effetti negativi di frane e altri rischi naturali. Le metodologie a disposizione sono di varia natura. In questo lavoro verranno presi in esame tre differenti modelli di valutazione della stabilità di tipo deterministico utilizzando i dati geologici provenienti dal progetto CARG (Cartografia Geologica d’Italia scala 1:50.000) di ISPRA (Istituto Superiore per la Protezione e la Ricerca Ambientale) per la realizzazione del foglio geologico dell’area di Antrodoco (RI) in fase di pubblicazione. I risultati della modellazione verranno valutati attraverso il calcolo di indici sintetici e la comparazione di curve ROC (Receiver Operating Characteristic) relative ai diversi modelli implementati.
In questo lavoro si intende mettere in risalto la variazione del fattore di sicurezza in funzione sia dell’utilizzo di differenti modelli per la verifica della stabilità che di differenti modelli per la valutazione dell’infiltrazione.

Stampa pdf rivista

Introduzione

La valutazione della suscettibilità da frana di un territorio può essere realizzata utilizzando un approccio di tipo deterministico ovvero attraverso metodologie che prevedono la modellazione fisica del fenomeno. Tali modelli sono caratterizzati da diversi livelli di complessità potendo essere monodimensionali, bidimensionali o tridimensionali e utilizzare sistemi di equazioni più o meno complessi per la risoluzione delle leggi fisiche che descrivono le condizioni di rottura, spostamento e arresto dei corpi in frana (Montgomery and Dietrich, 1994; Iverson, 2000; Baum et al., 2008; Lu and Godt, 2008; Baum and Godt, 2010). Per la soluzione di queste leggi viene richiesta la conoscenza di parametri geotecnici o geomeccanici puntuali, specifici per le diverse tipologie di frana (crolli, scivolamenti, colate, etc.) ottenibili necessariamente attraverso un rilevamento in sito e prove di laboratorio. I risultati della modellazione sono fortemente condizionati dalla qualità e dal livello di dettaglio di tali parametri di input. Per queste ragioni, l’analisi deterministica viene applicata su siti specifici e per fenomeni particolari (scala di versante), e generalmente non è molto utilizzata per effettuare una valutazione spaziale della pericolosità a piccola scala (Fell et al, 2008).
I modelli per l’analisi di suscettività da frana sono solitamente composti da un modulo geo-meccanico (che nel caso più comune è riconducibile alla stabilità di un Pendio Infinito) e uno idrologico accoppiati.
Il modello idrologico di infiltrazione è un aspetto fondamentale per la valutazione della stabilità e deve essere in grado di simulare processi quali l’infiltrazione, lo scorrimento superficiale, il flusso sub-superficiale e convertire l’altezza di pioggia in altezza del fronte saturo.
I modelli attualmente in uso possono differenziarsi per l’aspetto idrologico oppure per il fatto di essere stazionari o transitori.
Nel presente lavoro sono stati applicati tre modelli per l’analisi di stabilità distribuita:
a)    modello stazionario per un pendio infinito (Hammond et al, 1992);
b)    modello dell’afflusso critico (Pack et al, 1999);
c)    modello transitorio che utilizza la formulazione di Iverson (2000).
L’utilizzo di differenti metodi di analisi su una stessa area di studio e il confronto tra i risultati ottenuti può contribuire a migliorare la qualità e l’affidabilità di ciascun modello, a identificare i principali fattori predisponenti e innescanti dei fenomeni analizzati, a trascurare gli aspetti meno determinanti e, infine, a scegliere la metodologia più appropriata per raggiungere un obiettivo specifico. (Zizioli et al, 2013)
Tutte le valutazioni sono state effettuate utilizzando i sistemi informativi territoriali che, mediante operazioni di overlay, permettono di costruire le relazioni algebriche su più livelli informativi costituenti i parametri di input dei modelli considerati. L’unità territoriale di base utilizzata per l’elaborazione è un pixel avente dimensioni 20×20 m.

Area di studio

L’area di studio del presente lavoro è quella compresa nell’area di Antrodoco della Carta Geologica d’Italia scala 1:50.000 (Progetto CARG), ha una superficie di circa 650 km2 suddivisa tra le province di Rieti e L’Aquila, il territorio è prettamente montano e comprende le cime del monte Terminillo, che supera i 2000 metri s.l.m. Il territorio ricade nell’Appennino centrale ed il suo particolare interesse geologico è legato alla presenza di successioni sedimentarie che abbracciano paleoambienti che variano dal bacino alla piattaforma interna, con tutte le facies di passaggio, con età dal Triassico superiore al Neogene (Pantaloni et al. 2004; Berti et al. 2009).
Per definire gli aspetti litologici è stata utilizzato il Foglio n. 139 “L’Aquila” della Carta Geologica d’Italia alla scala 1:100.000 (SGI, 1955), che rimane al momento l’unica fonte geologica ufficiale del territorio nazionale, in attesa del completamento della cartografia geologica alla scala 1:50.000. La qualità del dato risulta piuttosto bassa, a causa della scala di rilevamento e del mancato aggiornamento della cartografia, che risale ai primi anni ’50.
Sono stati effettuati degli accorpamenti su litologie con analogo comportamento meccanico, in particolare su formazioni carbonatiche a spiccato comportamento litoide, definendo in tal modo 8 classi (Depositi quaternari, Depositi quaternari cementati, Flysch, Calcari marnosi, Marne e Marne argillose, Alternanze di calcari, argille e marne, Calcari stratificati e Dolomie) la cui distribuzione è riportata in fig. 1.

 

Fig. 1 - Mappa litologica e distribuzione fenomeni franosi: (1) Depositi quaternari; (2) Depositi quaternari cementati; (3) Flysch; (4) Calcari marnosi; (5) Marne e Marne argillose; (6) Alternanze di calcari Argille e marne; (7) Calcari Stratificati; (8) Calcari dolomitici  e Dolomie; (9) Inventario frane.

Per quanto riguarda l’archivio dei dissesti sono stati integrati i dati relativi all’Inventario dei Fenomeni Franosi Italiani (IFFI) (ISPRA, 2008) attraverso verifiche e osservazioni di campagna, rilievi foto interpretativi e indagini storiche su archivi.
Da tale analisi è stata realizzata la carta della distribuzione dei dissesti, con un grado di dettaglio alla scala 1:10.000.
In totale sono state censite 366 aree in frana comprendenti anche elementi codificati come scarpate ad alto angolo soggette a crollo e Deformazioni Gravitative Profonde di Versante (DGPV). Questi elementi areali sono stati esclusi dall’elaborazione ottenendo quindi un dataset di 214 frane per una superficie complessiva di circa 4,5 km2. Le frane interessano prevalentemente le successioni carbonatiche di piattaforma, scarpata e bacino (58%) e la formazione del flysch della Laga (29%). In misura minore interessano le successioni marnoso-arenacee del Miocene (5%) e i depositi di detrito (8%).
Per quanto riguarda la tipologia del movimento, se si considera il numero degli eventi franosi, prevalgono le frane da crollo (56%), cui seguono gli scivolamenti traslativi (18%). Se si considera l’estensione delle aree in frana, ponendo a parte le DGPV (Deformazioni Gravitative Profonde di Versante), i fenomeni di crollo continuano a essere tra i più rappresentati, anche se superati in estensione dalle frane complesse.


Metodologie

Modello stazionario per un pendio infinito

Il più semplice modello utilizzato è di tipo stazionario, accoppia un modello idrologico come proposto da Beven e Kirby (1979) e uno di stabilità dei pendii all’equilibrio limite.
Le principali assunzioni alla base della formulazione sono:
pendenza costante del pendio;
superficie di rottura piana e parallela al pendio tipicamente localizzata all’interfaccia tra due litologie con resistenza molto differente;
criterio di resistenza del terreno alla Mohr-Coulomb;
flusso stazionario parallelo al pendio;
assenza di flusso tra substrato e coltre.
In particolare al crescere dell’altezza hw dello strato di terreno saturo corrisponde una diminuzione della tensione agente normalmente al piano di rottura e, di conseguenza, della resistenza al taglio disponibile.
Questo approccio fa riferimento alla stabilità di un pendio indefinito applicabile a frane di ridotto spessore in relazione alla estensione longitudinale (scivolamenti traslativi, colate, debris flows) trascurando il contributo di resistenza dato dalle parti laterali, sommitali e basali della coltre.
I processi di instabilità investigati interessano spesso i terreni sciolti della copertura eluvio-colluviale, contraddistinti da una associazione con le acque meteoriche; tali dissesti si manifestano come uno scivolamento di suolo ed evolvono in colate entro le incisioni torrentizie.
È evidente che fenomeni roto-traslativi non possono essere adeguatamente modellati con questo algoritmo.
Il fattore di sicurezza per un pendio infinito per un materiale dotato di coesione ed in presenza di falda è il seguente:

dove:
c = coesione
γ = peso di volume del terreno
γw = peso di volume dell’acqua
h = profondità della superficie di rottura
hw = altezza piezometrica in corrispondenza alla superficie di rottura
β = pendenza del versante
φ = angolo di attrito sulla superficie di rottura

In questa relazione i parametri che sono necessari sono di tre tipi differenti:
    caratteristiche meccaniche, c, φ, γ ottenuti da una caratterizzazione geomeccanica;
    caratteristiche geometriche, β dal modello digitale del terreno e h spessore massimo dal rilevamento geologico;
    hw altezza piezometrica che è ottenuta dal modello idrologico utilizzato ed è funzione delle caratteristiche idrauliche e dall’apporto meteorico.
I valori utilizzati sono sintetizzati nella tabella 1.

Tab 1 - Parametri geotecnici dei litotipi utilizzati nei modelli.

Il modello idrologico di Beven & Kirby (1979) prevede un meccanismo di saturazione che procede dal basso verso l’alto e la definizione di un’area di drenaggio per unità di larghezza misurata lungo le isoipse.
Questo valore è definito dal modello idrologico secondo le seguenti assunzioni:
il deflusso segue i gradienti topografici ricavati dal modello digitale del terreno e dipende dalla relativa area di contributo a;
il deflusso è in equilibrio con un apporto dato dalla piovosità efficace R;
la capacità di deflusso è pari a T sinθ, dove T è la trasmissività del suolo.
L’altezza relativa di acqua nella coltre potenzialmente instabile si definisce come:

 


 

Attraverso funzioni di map algebra del software ArcGis (Esri) e in particolare attraverso lo strumento “raster calculator”, è possibile determinare il fattore di sicurezza per pixel, sotto le ipotesi precedentemente definite.
I risultati ottenuti mostrano le aree a maggiore propensione all’instabilità in funzione del fattore di sicurezza, visibili nella mappa di figura 2, dove sono state identificate 4 classi a pericolosità decrescente (Molto alta, alta, media e bassa). I limiti adottati per definire le classi sono i seguenti: <1, 1.0-1.5, 1.5-3, >3.

Fig.2 - Mappa di suscettibilità ottenuta dal modello stazionario

Afflusso critico (Shalstab)

Tale modello (Dietrich & Montgomery 1998; Montgomery & Dietrich 1994) utilizza la medesima teoria di quello precedentemente esposto, ovvero di stabilità di un pendio infinito con una coltre che si mobilizza su un basamento stabile.
A differenza del precedente, l’output del modello è costituito dalle altezze critiche di afflusso meteorico che generano instabilità.
Combinando infatti le equazioni di Darcy con quella relativa alla valutazione della stabilità di un pendio infinito si perviene alla seguente espressione:
 


In cui è rilevante il rapporto q/T che ingloba caratteristiche idrauliche rappresentate dalla trasmissività T e l’apporto stazionario q inteso come afflusso meteorico critico tale da instabilizzare la cella; le altre variabili geometriche sono contemplate nel fattore a/bsinθ.
Nella figura 3 viene esplicitata la simbologia adottata.

 

Fig - 3 SHALSTAB modello concettuale secondo Montgomery and Dietrich, 1994).   a) pendio infinito; b) modello idrologico (p, precipitazione; e, evapotraspirazione; r, drenaggio profondo; a, area drenata; h, altezza del livello dell’acqua; z, spessore del suolo; u, velocità del flusso subsuperficiale; hw, livello dell’acqua del flusso superficiale; θ, pendenza del versante; b, larghezza del canale di flusso).                   Fonte: “Recommended procedure for validating landslide model 7th Framework Programme Safeland D2-8 2011” (Modificata)

Il livello di saturazione necessario per la valutazione della stabilità risulta dal prodotto di questi fattori; si può pertanto esprimere la stabilità in funzione dei predetti parametri idraulici e naturalmente di quelli meccanici.
Esplicitando il valore di q si ottiene:



Come indicato dagli autori è stato assunto un intervallo temporale pari a 24 ore necessario a raggiungere le condizioni di stazionarietà.
Il metodo per la determinazione delle direzioni di deflusso fa riferimento all’algoritmo Multiple Flow Direction di Quinn et al. (1991) che verrà successivamente definito.
Come nel caso precedente, utilizzando lo strumento di calcolo “raster calculator”, è stato possibile determinare il valore di piovosità critica per pixel che produce instabilità.
Nell’ipotesi che l’afflusso meteorico non subisca perdite nella sua trasformazione in afflusso efficace alla falda, q rappresenta la variabile predittiva del modello, chiamata “pioggia critica”, intesa come quantità di pioggia richiesta per attivare un processo di franamento superficiale, ovviamente i valori più bassi indicano una maggiore propensione all’instabilità mentre i valori più alti indicano una maggiore propensione alla stabilità.
La dimensione temporale dei processi non viene colta esplicitamente dalla metodologia.
I risultati ottenuti sono stati sovrapposti alle frane effettivamente rilevate al fine di valutare la modellazione.
Nella figura 4 vengono mostrati i risultati ottenuti in termine di infiltrazione critica espressa in mm/day, necessaria per produrre instabilità.

 

Fig 4 - Mappa di suscettibilità ottenuta dal modello dell’afflusso critico

Secondo quanto proposto dagli autori sono identificate delle zone stabili e instabili a priori (indipendentemente dalle infiltrazioni) in funzione della pendenza e delle caratteristiche meccaniche, le rimanenti possono essere stabili o instabili a seconda dell’infiltrazione.
Per potere realizzare una mappa di pericolosità analoga a quanto già prodotto con 4 classi da bassa ad alta, non avendo come output dalla modellazione un valore di Fattore di Sicurezza, si è fatto uso del  metodo delle curve ordinatrici proposto da Chung e Fabbri (2003).
La definizione delle classi è stata effettuata in funzione del valore di Ratio of Effectiveness, che per essere significativo ed efficace ai fini predittivi deve avere valori molto grandi oppure prossimi allo zero. Secondo gli autori Chung e Fabbri (2003),  questo valore deve essere almeno maggiore di 3 o inferiore al massimo a 0,2, ma  preferibilmente maggiore di 6 o inferiore a 0,1.
Sono state pertanto individuate 4 classi di pericolosità aventi i seguenti limiti: <5 mm/day, 5-100 mm/day, 100-300 mm/day e >300 mm/day, corrispondenti a pericolosità molto alta, alta, media e bassa.

Modello transitorio (TRIGRS)

Un’ulteriore evoluzione di questi modelli è rappresentato da TRIGRS (Transient Rainfall Infiltration and Grid based Regional slope stability model; Baum et al 2008).
Tale modello utilizza un software proprietario, sviluppato in Fortran, che permette l’analisi del cambiamento della pressione dei pori e della variazione del fattore di sicurezza nel tempo dovuto all’infiltrazione meteorica, usando il metodo proposto da Iverson (2000).

 

Fig 5 - TRIGRS modello concettuale. Fonte: “Recommended procedure for validating landslide model 7th Framework Programme Safeland D2-8_2011”.

Il programma utilizza i seguenti parametri di input che vengono discretizzati cella per cella: intensità della precipitazione, pendenza, spessore della coltre, altezza dell’acqua iniziale, conducibilità idraulica verticale, diffusività idraulica, coesione, angolo di resistenza al taglio e peso di volume del suolo. Il software combina modelli di infiltrazione, di flusso subsuperficiale e ruscellamento e di stabilità per valutare gli effetti di eventi meteorici critici su vaste aree. La procedura utilizza un modello di stabilità e di infiltrazione come successivamente specificati.
a) Modelli di infiltrazione
Il processo è modellato attraverso la sovrapposizione di una componente transitoria e una stazionaria che dipende dalle condizioni iniziali.
Il modello di infiltrazione è basato sulla soluzione linearizzata di Iverson della equazione di Richards (1931).
La soluzione adottata per ottenere il carico idraulico al tempo t e alla profondità Z è:

 

Dove
Ψ = carico relativo alla pressione dell’acqua
t = tempo
Z = z/cosδ con Z profondità e z profondità normale all’angolo di pendenza δ
d = profondità dell’acqua misurata nella direzione z
IZLT = flusso superficiale stazionario
β = cos2δ - [IZLT /KZ]
Kz = conduttività idraulica nella direzione Z
Iz = flusso stazionario iniziale
InZ = flusso superficiale di data intensità per l’intervallo n
D1 = D0cos2δ dove D0 è la diffusività idraulica
N = numero di intervalli
H(t-tn) = funzione di Heavyside

Il primo termine dell’equazione rappresenta la parte stazionaria il rimanente quella transitoria.

b) Modello di flusso e ruscellamento
L’infiltrazione I quindi è la somma della Precipitazione P più un eventuale ruscellamento proveniente da monte Ru con la condizione di non superare la conduttività idraulica Ks secondo la seguente relazione:
I= P+ Ru  se  P+ Ru≤ Ks

Da ogni cella dove P+Ru è maggiore di Ks si genera un ruscellamento Rd che viene indirizzato alle celle sottostanti seguendo le seguenti condizioni:
Rd= P+ Ru− Ks se P+Ru-Ks ≥0
o
Rd =0 se P+Ru-Ks<0

c) Modello di stabilità
Come nei casi precedenti viene utilizzata la relazione per la valutazione della stabilità di un pendio infinito



dove     
c = coesione
φ = angolo di resistenza
γw = peso di volume dell’acqua
γs = peso di volume del terreno
β = pendenza
Z = spessore della coltre

L’ elaborazione prevede questi step:
inizializzazione del file contenente il DTM attraverso un modulo del programma denominato TopoIndex necessario per il funzionamento della routine di runoff in TRIGRS. Questa routine genera una griglia di fattori di peso che determinano quale proporzione del ruscellamento è trasferita a ciascuna cella adiacente. Dopo una serie di tentativi il sistema che ha dato migliori risultati tra tutti quelli possibili è quello del DINF di Tarboton (1997) che assume che il flusso proceda dalla pendenza maggiore, calcola la direzione del pendio e ripartisce il flusso attraverso due celle con la pendenza maggiore.

Definizione dei dati di input.

Le caratteristiche geotecniche ed idrauliche delle coltri sono le stesse utilizzate nelle precedenti simulazioni. Nella tabella 2 vengono sintetizzati i parametri di Permeabilità (K) e Diffusività (D0), non considerati dai precedenti modelli.

Tab.2 - Valori di Permeabilità (K) e Diffusività (D0) utilizzati nel modello.

Per quanto riguarda l’altezza iniziale dell’acqua è stata assunta coincidente con il bedrock, ipotizzando che precedentemente all’evento le coltri allo stato di pre-evento fossero asciutte.
In assenza di informazioni di dettaglio i parametri idraulici necessari alla modellazione sono stati determinati attraverso una analisi parametrica avente per finalità la massimizzazione del rapporto celle instabili correttamente classificate/celle instabili totali. Il fattore di sicurezza viene calcolato in diversi intervalli temporali e per diverse profondità della coltre. Tra tutti gli output possibili, relativi a diversi step temporali e per differenti spessori, nella figura 6 è riportato il Fattore di sicurezza minimo.
Sono state effettuate diverse elaborazioni tramite le quali è stato possibile valutare l’evoluzione dell’evento, nell’analisi dei risultati si fa riferimento a quella finale delle 24 ore.
I limiti delle classi di pericolosità sono gli stessi utilizzati nel modello del pendio infinito.

Confronto tra modelli di stabilità

La valutazione dell’affidabilità e della performance di un modello costituisce un aspetto fondamentale nella modellazione della pericolosità per frana.
Il criterio più rilevante per la valutazione della qualità è la valutazione della precisione del modello, che viene eseguita analizzando l’accordo tra i risultati del modello e dei dati osservati (Frattini et al 2010, Corominas & Mavrouli, 2011).
Tale confronto è realizzato attraverso la costruzione della matrice di confusione che, confrontando i valori reali (le frane censite nell’inventario) e quelli prodotti dal modello, permette il conteggio dei veri negativi (TN), veri positivi (TP), falsi negativi (FN) e falsi positivi (FP). L’unità spaziale utilizzata è il pixel pertanto il lavoro sopraccitato consiste nel conteggio dei pixel correttamente e non correttamente classificati. Un pixel valutato instabile dal modello, se sovrapposto ad una frana censita, dà luogo ad un TP, altrimenti ad un FP. Viceversa, un pixel valutato stabile, se sovrapposto ad una frana censita, dà luogo ad un FN, altrimenti ad un TN. Molti sono gli indici statistici sintetici che vengono utilizzati, in questo lavoro si intende presentare i risultati di due indici che utilizzano tutti i 4 gli elementi della tabella di contingenza. In particolare è stato calcolato l’Odd ratio skill score (Yule’s Q; Yule, 1900) attraverso la relazione:




tale indice è collegato all’Odds Ratio (Stephenson, 2000), che misura il rapporto di Odds tra predizioni vere e false, ma normalizzato ed oscillante tra i valori di +1 e -1.
Inoltre è stato calcolato l’Heidke’s skill score (k di Cohen; Heidke, 1926) che misura il grado di accordo tra due valutatori:



dove  Pobs è la proporzione dell’accordo osservato (ottenuta sommando i valori della diagonale principale della matrice confusione e dividendo per il numero totale dei casi), mentre Pexp è la proporzione dell’accordo dovuto al caso (ricavata sommando i totali marginali di riga e di colonna diviso il numero totale dei casi).
Per costruire la matrice confusione è necessario inserire delle soglie che discriminano i pixel stabili da quelli instabili. Nel caso dei modelli di stabilità che forniscono come risultato un fattore di sicurezza, tale soglia è pari ad 1, mentre nel caso del modello di afflusso critico sono state costruite due matrici confusione utilizzando due diversi valori di cut-off:
 a) 5 mm/day, soglia che si ottiene dalla curva ordinatrice tale che il Ratio of Effectiveness sia significativo (Chung e Fabbri 2003);
 b) 100 mm/day, desunta dalla curva ROC (descritta successivamente), considerando il punto più prossimo allo spigolo in alto a sinistra.
Dalla sola analisi di questi valori si dovrebbe osservare che il migliore metodo risulta essere quello Transitorio, con un valore di Q pari a 0.7 piuttosto soddisfacente.

 

MODELLO Yule’s Q

Cohen’s K

STAZIONARIO (FS=1)

0,56 0,020

TRANSITORIO (FS=1)

0,70 0,021

AFFLUSSO CRITICO (AC=5)

0,42 0,016
AFFLUSSO CRITICO (AC=100) 0,36 0,006

 

Tab. 3 - Valori di k di Cohen e del Q di Yule calcolati sui modelli di stabilità; tra parentesi il valore di soglia considerato per generare la matrice confusione.

In realtà dall’analisi dei valori di k di Cohen è evidente un’altra situazione: le differenze tra i diversi tipi di modellazione sono molto più sfumate e soprattutto i bassi valori che si evidenziano sono sintomo di una bassa performance in termini di concordanza tra dato simulato e reale. Questa osservazione trova riscontro quando, anziché utilizzare un indicatore che ha la limitazione di considerare un solo valore di cut-off, si stima la performance attraverso la definizione di una curva ROC (Figura 7) utilizzando differenti valori di cut-off (Begueria, 2006) e calcolando il parametro AUC (Area Under Curve) come misura della qualità complessiva del modello (Hanley and McNeil, 1982).

Fig. 7 - Curve ROC relative ai tre modelli di stabilità (curva A - modello stazionario; curva B - modello transitorio; curva C - modello dell’afflusso critico)

In modo controintuitivo, la modellazione che utilizza Shalstab, avente i più bassi valori di  k e Q, è caratterizzata da valori di AUC migliori degli altri modelli (AUC=0.63, curva C di figura 7); il modello stazionario e quello transitorio hanno invece valori simili ma leggermente inferiori, rispettivamente 0.60 (curva A) e 0.59 (curva B).
Dalla forma delle curve si evince che gli ultimi due modelli sono molto simili, risultato in linea con il fatto che utilizzano lo stesso modello di stabilità e pertanto in condizioni di totale saturazione anche se ottenuta da modelli di infiltrazione differenti tendono a dare risultati analoghi. È evidente che dati i valori di AUC la performance di queste modellazioni è medio-bassa, risultato che è in concordanza con i bassi valori del k di Kohen.
È evidente che i modelli proposti non corrispondono in modo adeguato alle tipologie di fenomeni avvenuti ma sicuramente più rilevante ancora è la bassa risoluzione del DTM, la caratterizzazione geotecnica insufficiente per un’areale così vasto e soprattutto la definizione delle regioni basate su una carta litologica avente una scala non di dettaglio.
I risultati dell’analisi suggeriscono che, nell’area di studio per la determinazione della suscettività alle frane, il contributo di un modello fisicamente basato, completo di componente geo-meccanica e componente idrologica accoppiate, è assolutamente trascurabile rispetto ad un semplice modello geo-meccanico basato sulla pendenza e sulle caratteristiche meccaniche come quello del Pendio Infinito. La variabile pendenza contiene quindi quasi tutta l’informazione disponibile, come confermato dall’implementazione di alcuni modelli statistici sull’area di studio.
I modelli che ipotizzano, nella componente idrologica, un movimento dell’acqua sub-superficiale prevalentemente parallelo al versante come SHALSTAB non sono adatti alle permeabilità molto basse caratteristiche delle argille. I modelli che ipotizzano il flusso verticale come TRIGRS, basato sul lavoro di Iverson che si basano su assunzioni più adatte a litologie a bassa permeabilità, dovrebbero dare migliori risultati, ma come si è mostrato il contributo aggiuntivo che possono dare allo studio di suscettività (quindi solo alla localizzazione spaziale) è trascurabile rispetto al solo Pendio Infinito.

Modelli di infiltrazione

È fondamentale nella valutazione del fattore di sicurezza in una coltre potenzialmente instabile parzialmente satura valutare l’altezza relativa di acqua nella coltre stessa.
A parità di modello di calcolo utilizzato gli algoritmi con cui vengono tramutati gli afflussi meteorici in altezza d’acqua che satura la coltre possono essere differenti. Come già espresso nella (2) che utilizza il semplice modello di Beven & Kirby dipende da aspetti geometrici ed idraulici, considerando che i valori di pendenza e trasmissività possono essere considerati costanti, l’unico fattore che può essere oggetto di variazione è l’area contribuente.   
L’area contribuente è la porzione planimetrica di area a monte di un elemento topografico discreto, estesa sino allo spartiacque, che contribuisce alla formazione del deflusso nell’elemento considerato. Il calcolo dell’area contribuente non è univoco, ma dipende dall’adozione di alcuni algoritmi per l’identificazione delle direzioni di deflusso. In questo lavoro sono stati utilizzati quattro differenti algoritmi, implementati in software GIS, a partire dal più antico ai più recenti e sofisticati; le quattro distribuzioni di area contribuente ricavate sono state utilizzate per la determinazione del fattore di sicurezza distribuito considerando un modello di stabilità che fa riferimento al pendio infinito.
Le elaborazioni sono state effettuate utilizzando il software GIS SAGA (System for Automated Geoscientific Analyses; Conrad et al, 2015; Olaya & Conrad, 2008), in particolare la libreria Hydrology con il modulo Catchment Area (Top-Down). Il DEM utilizzato è stato preprocessato eliminando le depressioni esistenti.
Di seguito vengono esposti i quattro algoritmi utilizzati:
a) Algoritmo D8
Nell’algoritmo D8 (O’Callaghan e Mark, 1984) viene ipotizzato che il deflusso che attraversa una generica cella del DEM (modello digitale delle elevazioni) si muova interamente lungo la direzione che esprime il massimo gradiente topografico negativo tra la cella in esame e quelle adiacenti. Le possibili direzioni esplorate dall’algoritmo sono otto ovvero quelle degli assi di simmetria che si irradiano dal centro della cella (da cui il nome dell’algoritmo D8).
b) Algoritmo Dinf
Secondo l’algoritmo D-infinito (Tarboton, 1997) il deflusso che attraversa una certa cella si muove verso valle lungo una sola direzione, ma a differenza di quanto accade per il D8, tale direzione non deve necessariamente coincidere con uno degli otto assi principali di simmetria. La superficie di ciascuna cella viene suddivisa in otto porzioni triangolari con vertice condiviso al centro della cella (figura 8a).
 

 


Fig 8 - Flow direction secondo il modello D∞: (a) suddivisione del flusso in due celle adiacenti (b) calcolo della pendenza nella singola faccetta. Fonte: “Comparison of Different Grid Cell Ordering Approaches in a Simplified Inundation Model_Yang et al_2015” (modificata)

a ciascuna di esse viene assegnata una pendenza in funzione della quota dei baricentri delle celle confinanti. Con riferimento a una sola delle porzioni triangolari, la pendenza è rappresentata dalle seguenti formulazioni:



in cui ei sono le quote delle celle e di le distanze (figura 8b).
Noto il vettore (s1, s2), il modulo s della pendenza e la direzione r associata alla porzione di superficie sono definite dalle relazioni:


Se la direzione r risulta esterna alla superficie della porzione triangolare rappresentata in figura, la direzione di deflusso viene fatta coincidere con il più prossimo dei due lati e0e2, e0 e1 ed il deflusso ridirezionato alla cella adiacente in direzione cardinale o diagonale senza subire ripartizione. Nel caso in cui la direzione individuata ricada all’interno della porzione, essa viene ripartita tra le due celle in direzione cardinale e diagonale in modo inversamente proporzionale agli angoli α1 ed α2 formati da r con le due direzioni principali di simmetria. Ne consegue che il deflusso che caratterizza una cella, nel muoversi verso le celle adiacenti a quota inferiore, o non viene ripartito o viene ripartito al più tra due di esse. La pendenza locale associata alla cella è la maggiore tra quelle calcolate per le otto porzioni triangolari.
c) Algoritmo Multiple Flow
L’algoritmo Multiple Flow (Quinn et al., 1991) è, tra quelli considerati, quello che presenta la maggior dispersione del deflusso tra una generica cella e quelle adiacenti a quota inferiore (fig. 9). In esso non viene effettuata la “scelta” di una direzione, ma si ipotizza che il deflusso si ripartisca tra tutte le celle sottostanti in misura proporzionale a: i) pendenza tra cella ricevente e cella contribuente, ii) lunghezza del tratto di isoipsa che il deflusso attraversa ortogonalmente per raggiungere la cella ricevente (Vedi figura sottostante).

 


 

Fig 9 - Flow direction using the multiple flow direction model. Fonte: ‘The prediction of hillslope flow paths for distributed hydrological modelling using digital terrain models’ -Quinn 1991

Come per il D8, le possibili direzioni di deflusso sono rappresentate dagli assi principali di simmetria della cella: il deflusso può essere ripartito al massimo in otto parti, se la cella risulta essere la più alta tra tutte le circostanti. Le pendenze tra la cella considerata e le adiacenti a quota inferiore fungono da pesi per determinare la porzione di deflusso in ciascuna delle direzioni. Se si associa valore unitario all’area contribuente relativa ad una singola cella, la frazione di deflusso ΔFacci che viene attribuita alla i-esima cella ad essa adiacente a quota inferiore è data da:



dove
    n = numero totale di celle adiacenti di quota inferiore;
    tan βi = pendenza verso l’i-esima cella;
    Li = lunghezza della linea di livello nella i-esima direzione.
La pendenza associata alla cella centrale risulta, in maniera analoga, dalla media delle pendenze verso le celle adiacenti a quota inferiore pesate dalle rispettive lunghezze dei tratti di isoipsa che intercettano il deflusso:




d) Algoritmo Maximum Downslope Gradient
Tale algoritmo, di seguito chiamato MDG, (Qin et alii, 2007) rappresenta l’evoluzione del metodo Multiple Flow implementato da Quinn, un approccio adattivo per determinare la ripartizione di flusso sulla base di attributi topografici locali. In questo approccio, la ripartizione flusso è modellata da una funzione di ripartizione che si basa sulla massima pendenza locale.
Gli autori ritengono che le condizioni del terreno locali controllano la ripartizione del flusso in ogni cella. È più ragionevole utilizzare condizioni topografiche locali per determinare la ripartizione del flusso che usare attributi topografici globali (come ad esempio il flow accumulation).
Utilizzando le condizioni topografiche locali l’algoritmo è più adattabile alle condizioni che controllano direttamente la ripartizione del flusso locale. Questa idea comporta lo sviluppo di un nuovo approccio alla ripartizione flusso secondo MFD. Ci sono due fasi nella costruzione di questo nuovo approccio.
La prima è la definizione della caratteristica topografica locale (e) che descrive direttamente l’effetto del terreno sulla ripartizione del flusso. La seconda fase è la costruzione di una funzione di e (f (e)) che serva a calcolare la ripartizione di flusso che si adatta alle condizioni del terreno locale.

Confronto tra modelli di infiltrazione

In virtù della maggiore “dispersione del deflusso” da una cella a quelle di valle che caratterizza l’algoritmo Multiple Flow, le distribuzioni di area drenata ottenute da esso risultano organizzate in modo più graduale che quelle derivanti dai precedenti algoritmi.
Nelle figure 10 - 13 vengono riportate quattro mappe raffiguranti l’area drenata espressa in forma logaritmica calcolata rispettivamente mediante l’algoritmo D8, DINF, MFD e MDG.

 

Fig 10 - Mappa dell’area contribuente ottenuta con algoritmo D8

 

Fig 11 - Mappa dell’area contribuente ottenuta con algoritmo DINF

 

Fig 12 - Mappa dell’area contribuente ottenuta con algoritmo MDG

Fig 13 - Mappa dell’area contribuente ottenuta con algoritmo MFD

Si nota una netta diversità nella distribuzione delle classi di area drenata (Figura 14),

 

Fig 14 - Distribuzione dell’area contribuente nei quattro modelli di infiltrazione utilizzati.

indice del fatto che le condizioni di calcolo sono effettivamente differenti, inoltre è evidente una maggiore connessione del reticolo idrografico negli ultimi due algoritmi. Utilizzando il modello di stabilità che fa riferimento al pendio infinito è stata effettuata la valutazione del Fattore di sicurezza per tutti modelli di infiltrazione adottati. Al fine di valutare la performance di questi sono state determinate le curve ROC di cui si è detto in precedenza, visibili nella figura 15, dove per confronto è stato inserito anche il modello transitorio; i valori di AUC sono sintetizzati nella tabella 4. In termine di Fs le mappe non differiscono in maniera rilevante indizio che il meccanismo di infiltrazione non gioca un ruolo rilevante nell’area di studio considerata.

MODELLO AUC
D8 0.6033
DINF 0.6101
MFD 0.6134
MDG 0.6101
TRIGRS 0.5923

 

Tab. 4 - Valori di AUC dei modelli di infiltrazione considerati

 

Fig. 15 - Curve ROC dei modelli di infiltrazione considerati

Conclusioni

Sono stati messi a confronto 3 tra i più utilizzati metodi deterministici per la valutazione della stabilità distribuita valutando le performance attraverso indici sintetici e curve ROC. Nel lavoro si è dimostrato che in prima analisi il modello che ha presentato migliori risultati sembra essere quello transitorio (TRIGRS) con valori di Yule’s Q pari a 0.7, nettamente superiore a quello degli altri due metodi; dal valore assunto dal k di Cohen, tuttavia, è risultato evidente che la corrispondenza tra simulato e reale è piuttosto bassa. La bassa performance è dimostrata dal valore di AUC delle curve ROC che raggiungevano nel caso migliore il valore di 0.6. Il test effettuato utilizzando quattro differenti modelli di infiltrazione ha permesso di osservare che nonostante la definizione dell’area di contributo abbia dato risultati anche piuttosto diversi, la valutazione in termini di Fattore di sicurezza non presenta significative differenze.
Le curve ROC effettuate infatti presentano andamenti molto simili tra loro come anche gli AUC che si attestano sempre su valori piuttosto modesti. Alla luce dei risultati ottenuti si pone l’attenzione sull’utilizzo dei modelli per l’analisi di suscettività e ancora prima in fase di scelta degli stessi: un modello che si basi sulla concentrazione topografica del flusso sub-superficiale potrebbe infatti condurre a interpretazioni fuorvianti se non attentamente calibrato e confrontato con i risultati raggiungibili senza considerare gli effetti di questa ipotesi.
Ovviamente i risultati del presente studio non escludono l’utilità della modellazione idrologica dell’intero versante: in tutti i casi in cui è necessaria un’analisi dei tempi di attivazione e delle soglie innescanti, non si può certo fare a meno dei modelli idrologici, soprattutto per eventi che interessano la coltre superficiale, dove la pioggia, e il conseguente incremento delle pressioni neutre, rappresenta il fattore scatenante. Esistono inoltre terreni con caratteristiche diverse dall’area esaminata nel presente lavoro, dove il flusso da monte non può essere trascurato e dove quindi la modellistica idrologica è necessaria anche per l’analisi di suscettività. Come citato in premessa, i risultati della modellazione sono fortemente condizionati dalla qualità e dal livello di dettaglio dei parametri di input; per tale motivo la loro applicazione risulta molto più appropriata in situazioni locali ove siano disponibili dati affidabili di tipo geotecnico, idrologico e si disponga di un DTM di dettaglio. Si intende sottolineare infine che in fase di validazione dei risultati è preferibile utilizzare metodi che utilizzano tutti i valori della matrice di confusione anziché utilizzare indici sintetici che possono portare a valutazioni non realistiche. In particolare le curve ROC e la relativa AUC forniscono una stima complessiva e sintetica della bontà della simulazione. Si sottolinea anche l’aspetto della difficoltà della stima della performance di un modello in cui la popolazione dei casi negativi (in questo caso il territorio stabile) risulta di gran lunga preponderante rispetto ai positivi. La concordanza tra il modello previsionale e la realtà infatti può essere anche dovuta al caso ed utilizzando metodologie che non ne tengano conto può determinare una affidabilità della simulazione non realistica.
Inoltre, a prescindere dal modello di calcolo adottato, dovendo rappresentare i risultati ottenuti con la consueta suddivisione in 4 classi esiste un’altra fonte di incertezza rappresentata dalla scelta dei limiti tra classi che vengono utilizzati, per i quali è consigliabile adottare una suddivisione il più possibile standardizzata come proposto p. es. da Chung e Fabbri (2003) mediante l’utilizzo del Ratio of effectiveness.

Bibliografia

Baum R.L., Savage W.Z., Godt J.W. (2008). TRIGRS – A FORTRAN program for transient rainfall infiltration and grid-based regional slope stability analysis, vers. 2.0, U.S. Geol. Survey Open-File Rep. 2008-1159, 75pp.
Begueria S. (2006). “Validation and evaluation of predictive models in hazard assessment and risk management”, Nat. Hazards, 37, 315– 329.
Beven, K. J. and Kirby, M. J.: A physically based variable contributing area model of basin hydrology, Hydrol. Sci. Bull, 24, 43–69, 1979.
CHANG-JO F. CHUNG  and ANDREA G. FABBRI  (2003) Validation of Spatial Prediction Models for
Landslide Hazard Mapping  Natural Hazards 30: 451–472,
Dietrich, W. E., and D. R. Montgomery, 1998, Hillslopes. channels and landscape scale. G. Sposito, editor, Scale dependence and scale invariance in hydrology, Cambridge University Press, Cambridge, England.
Montgomery, D. R. and Dietrich, W. E.: A physically based model for the topographic control of shallow landsliding,Water Resour. Res., 30, 1153–1171, 1994
Iverson, R.M., 2000, Landslide triggering by rain infiltration: Water Resources Research, v. 36, no. 7, p. 1,897–1,910.
Paolo Frattini, Giovanni Crosta , Alberto Carrara  (2010) Techniques for evaluating the performance of landslide susceptibility models.  Engineering Geology 111 62–72
Heidke, P. (1926). Berechnung der erfolges und der gute der windstarkevorhersagen im sturmwarnungdienst. Geogr. Ann., 8, 301- 349.
ISPRA Institute for Environmental Protection and Research (2008) Landslides in Italy, Special Report 2008. 83, ISPRA, Italy
Montgomery D.R., Dietrich W.E. (1994). “A physically based model for the topographic control of shallow landsliding”, Water Resour. Res., 30, 1153–1171.
Pack R.T., Tarboton D.G., Goodwin C.G. (1999). “SINMAP 2.0 – A Stability Index Approach to Terrain Stability Hazard Mapping, User’s Manual”. Forest Renewal B.C. under Research Contract No: PA97537-0RE.
Vercesi P., Scagni G. (1984). “Osservazioni sui depositi conglomeratici dello sperone collinare di Stradella”, Rend. Soc. Geol. It., 7, 23–26.
Tarboton, DG (1997) A new method for the determination of flow directions and upslope areas in grid digital elevation models. Water Resour Res 33: pp. 309-319
L A Richards  (1931) Capillary conduction of liquids through porous media Physics, 1 pp. 318–333
Zizioli D., Meisina C., Valentino R., Montrasio L. (2013). “Comparison between different approaches to modeling shallow landslide susceptibility: a case history in Oltrepo Pavese, Northern Italy”, Nat. Hazards Earth Syst. Sci., 13, 559–573.
Yule, G.U., 1900. On the association of attributes in statistics. Philosophical Transactions of the Royal Society of London 194A, 257–319.
Robin Fell, Jordi Corominas, Christophe Bonnard, Leonardo Cascini, Eric Leroi, William Z. Savage  (2008) “Guidelines for landslide susceptibility, hazard and risk zoning for land use planning” pubblicato su Engineering Geology, 102:85-98 – doi:10.1016/j.enggeo.2008.03.022.  
Corominas J. & Mavrouili O., 2011. Recommended Procedures for Validating Landslide Hazard and Risk Models and Maps. SafeLand FP7, Deliverable D2.8, 162 pp
Stephenson, D.B., 2000. Use of the “odds ratio” for diagnosing forecast skill. Weather Forecasting 15, 221–232.
Hanley, J.A., McNeil, B.J., 1982. The meaning and use of the area under a receiver operating characteristic (ROC) curve. Radiology 143 (1), 29–36.
Conrad, O., Bechtel, B., Bock, M., Dietrich, H., Fischer, E., Gerlitz, L., Wehberg, J., Wichmann, V., and Böhner, J. (2015): System for Automated Geoscientific Analyses (SAGA) v. 2.1.4, Geosci. Model Dev., 8, 1991-2007, doi:10.5194/gmd-8-1991-2015
Olaya, V., Conrad, O., 2008. Geomorphometry in SAGA. In: Hengl, T. and Reuter, H.I. (Eds), Geomorphometry: Concepts, Software, Applications. Developments in Soil Science, vol. 33, Elsevier, 293-308 pp.
O’Callaghan, J.F. Mark, D.M. (1984): ‘The extraction of drainage networks from digital elevation data’, Computer Vision, Graphics and Image Processing, 28:323-344
Quinn, P.F., Beven, K.J., Chevallier, P., Planchon, O. (1991): ‘The prediction of hillslope flow paths for distributed hydrological modelling using digital terrain models’, Hydrological Processes, 5:59-79
Qin, C. Z., Zhu, A. X., Pei, T., Li, B. L., Scholten, T., Behrens, T. & Zhou, C. H. (2011): ‘An approach to computing topographic wetness index based on maximum downslope gradient’, Precision Agriculture, 12(1), 32-43.
Pantaloni M, Capotorti F, D’Ambrogi C, Di Stefano R (2004) Geological guide of the 348 Antrodoco sheet ISPRA. Internal document, ISPRA, Italy.
Berti D et al (2009) La geologia del Foglio n 348 Antrodoco. Memorie Descrittive Carta Geologica d’Italia LXXXVIII: 134
SGI Servizio Geologico d’Italia (1955) Carta Geologica d’Italia Scala 1:100 000, Foglio n 139 L’Aquila
 

« Indice del n. 73