INDICE

1. INTRODUZIONE
2. DESCRIZIONE DEL MODELLO

3. APPLICAZIONE DEL MODELLO ALLA PIANA COSTIERA FRA IL TORRENTE TRIPESCE ED IL FOSSO DELLA MADONNA

4. CONCLUSIONI

Bibliografia

ELENCO FIGURE

Fig. 1 Distribuzione dei pozzi utilizzati per l’elaborazione del modello idrogeologico della Piana costiera
Fig. 2 Ubicazione delle sezioni elaborate per la costruzione del modello idrogeologico
Fig. 3 Interpretazione stratigrafica dei dati delle perforazioni (Alcune sezioni rappresentative: Sez2,Sez4,Sez5)
Fig. 4 Caratterizzazione idraulica dei litotipi e sezioni idrogeologiche.
Fig. 5a,b Discretizzazione dell’area modellizzata: (a) griglia degli elementi finiti in pianta, (b) in una sezione.
Fig. 5c Discretizzazione dell’area modellizzata: distribuzione dei parametri idraulici sul livello 10 (rappresentativo di uno degli acquiferi principali)
Fig. 6 Piezometrica Aprile 1996
Fig. 7 Piezometrica Settembre 1996
Fig. 8 Carta della variazione piezometrica dall’Aprile 1996 al Settembre 1996
Fig. 9 Diagramma di calibrazione
Fig. 10 Simulazione della piezometrica per lo strato 10 alla fine del periodo 1 (Estivo)
Fig. 11 Direzioni di flusso simulate per l’area Mazzanta-Foce Fiume Cecina
Fig. 12 Carte della conducibilità elettrica specifica. Area Mazzanta - Foce del F. Cecina
Fig. 13 Simulazione delle velocità di flusso per lo strato 10, periodo estivo
Fig. 14a,b (a) Idrogrammi simulati nei due periodi (estivo ed invernale) di alcuni pozzi rappresentativi della zona Mazzanta-foce F. Cecina. (b) Idrogrammi simulati nei due periodi (estivo ed invernale) di alcuni pozzi appresentativi della zona a Sud del Fiume Cecina
Fig. 15 Bilancio a zone per l’anno 1996

 

  1. INTRODUZIONE

Un modello è un utensile disegnato per rappresentare una versione semplificata della realtà. Data una definizione così ampia è evidente che si utilizzano modelli nella vita di ogni giorno. Per esempio una carta stradale è un modo di rappresentare una complessa rete di strade in forma simbolica così che è possibile testare varie strade su carta piuttosto che sbagliare durante la guida. Una carta stradale può essere considerato una sorta di modello poiché è un modo di rappresentare la realtà in forma semplificata. In modo simile i modelli relativi alle acque sotterranee sono rappresentazioni della realtà e, se dovutamente costruiti, possono essere strumenti previsionali per la gestione della risorsa idrica sotterranea. Usando modelli è possibile colaudare "virtualmente" vari schemi di gestione e predire gli effetti delle azioni che si intende intraprendere. Naturalmente la validità delle previsioni dipenderà da quanto il modello approssima le condizioni reali. Dati di partenza affidabili sono indispensabili quando si utilizza un modello a scopo previsionale.
Modelli matematici sul flusso delle acque sotterranee sono stati usati fin dal 1800. Un modello matematico consiste del set di equazioni differenziali che descrivono il movimento delle acque sotterranee. L’attendibilità delle previsioni (simulazioni) dipende da quanto bene il modello approssima la realtà. Normalmente le assunzioni necessarie per risolvere un modello matematico in via analitica sono troppo restrittive (ad esempio le soluzioni analitiche generalmente approssimano in maniera attendibile solo un mezzo omogeneo ed isotropo) così che per risolvere casi reali è preferibile utilizzare le tecniche numeriche. Dal 1960, quando i computer ad alta velocità si sono diffusi, i modelli numerici sono diventati il metodo favorito per studiare il movimento delle acque sotterranee ed il trasporto di contaminanti.

Nell’esempio di seguito riportato è stato utilizzato "MODFLOW": tale modello fu messo a punto dall’USGS negli anni ’70 è stato poi reso più maneggevole ed impiegabile anche da parte di non informatici negli anni successivi. Attualmente la vecchia versione in fortran ed una delle prime versioni per DOS – W 3.1 sono scaricabili gratuitamente da INTERNET. Da sottolineare che le attuali versioni di "Visual Modflow" usano ancora come base il metodo numerico dell’originale versione fortran.
Di seguito si riportano i caratteri essenziali del modello ad elementi finiti MODFLOW. 

 

MODFLOW, è un modello a differenze finite, che simula il flusso di acque sotterranee nelle tre dimensioni spaziali (X,Y,Z); esso incorpora i modelli bi e tridimensionali descritti da Trescott (1975), Trescott e Larson (1976), Trescott, Pinder and Larson (1976) e ampiamente utilizzati dall’ USGS. La sua implementazione è stata condotta in modo da garantirne la facile modificabilità, la semplicità di uso e soprattutto la compatibilità con una grande varietà di computers.
Il modello offre la possibilità di simulare:

MODFLOW permette la simulazione in regime di flusso stazionario o in regime di flusso transiente con possibilità, in quest’ultimo caso, di variare per periodi definiti le caratteristiche di cui ai punti precedenti. In questa breve esposizione si prenderà in considerazione solo il regime transiente, utilizzato nella simulazione per la bassa Val di Cecina.

2.1 Note Matematiche

Il movimento tridimensionale dell’acqua sotterranea (a densità costante) attraverso un mezzo poroso può essere descritto dall’equazione a derivate parziali:

d/dx (Kxx dh/dx) + d/dy (Kyy dh/dy) + d/dz (Kzz dh/dz) - W = Ss dh/dt (1)

dove:
Kxx, Kyy e Kzz sono i valori della conducibilità idraulica lungo le direzioni x,y,z che sono state assunte parallele agli assi maggiori della conduttività idraulica (Lt-1);
h è il carico idraulico;
W è il flusso volumetrico di acqua per unità di volume dell’acquifero e rappresenta gli apporti o le perdite di acqua (t-1);
Ss è lo ‘’storage’’ specifico del materiale poroso (L-1);
t è il tempo.

L’equazione (1), insieme con le condizioni al contorno di carico idraulico e con le condizioni di carico idraulico iniziali, costituisce la rappresentazione matematica di un sistema di flusso delle acque sotterranee. Eccetto che per rari casi molto semplici, l’equazione (1) non può essere risolta per via analitica; sono stati perciò sviluppati vari metodi di calcolo che forniscono soluzioni più o meno approssimate. Un approccio alla soluzione è il metodo delle differenze finite in cui il sistema continuo descritto dall'equazione (1) è rimpiazzato da un set finito di punti discreti sia nello spazio che nel tempo, e le derivate parziali sono sostituite da termini calcolati dalle differenze di carico idraulico in questi punti. Dal processo sopra descritto si deducono sistemi di equazioni algebriche lineari di differenze; la loro soluzione provvede valori di carico idraulico h agli specifici punti del sistema. Questi valori rappresentano una approssimazione alla distribuzione , variabile con il tempo, del carico idraulico rispetto a quella che si otterrebbe dalla soluzione analitica dell’equazione (1).

2.2 Convenzioni della Discretizzazione

In base a quanto definito nel paragrafo precedente, in MODFLOW il ‘’sistema acquifero’’ viene discretizzato nelle tre direzioni spaziali; ne risulta una suddivisioni in righe e colonne sul piano X,Y ed in strati in direzione Z. Nel complesso il sistema in studio viene suddiviso in un numero finito di celle; all’interno di ciascuna cella esiste un punto chiamato nodo in corrispondenza del quale viene calcolato il carico idraulico.
Esistono due convenzioni per la definizione delle celle rispetto ai nodi:

Nodi al centro del blocco. I blocchi formati da set di linee parallele sono le celle ed i nodi sono al centro delle celle stesse;
Punti centrati. I nodi sono i punti di intersezione dei set di righe parallele e le celle sono tracciate intorno ai nodi a metà strada fra i nodi.

Nel metodo delle differenze finite tutte le variabili vengono discretizzate e quindi è richiesta anche una discretizzazione del dominio tempo; il tempo di interesse viene quindi suddiviso in intervalli discreti ‘’time step’’.

2.3 Equazione della differenza finita

Lo sviluppo dell’equazione di flusso con il metodo delle differenze finite deriva dall’equazione di continuità: la somma di tutti i flussi in ingresso ed in uscita dalle celle deve essere uguale alla variazione di storage nella cella. Assumendo costante la densità dell’acqua, l’equazione di continuità che esprime il bilancio del flusso per una cella è:

S Qi = Ss Dh DV (2)

dove:

Qi = la velocità di flusso in una cella (L3T-1)
Ss = ‘’Storage’’ specifico; volume di acqua che può essere immagazzinata per variazione unitaria di carico idraulico (L-1)
DV = Volume di una cella (L3)
Dh = Variazione di carico idraulico nell’intervallo di tempo Dt.

Il termine a destra dell’equazione (2) è equivalente al volume di acqua immagazzinata in un intervallo di tempo Dt per una data variazione di carico idraulico Dh. L’equazione (2) viene quindi definita in termini di ingresso e guadagno in ‘’storage’’: uscite e perdite vengono rappresentate definendo le uscite come flussi negativi e le perdite come guadagni negativi. Flussi esterni all’acquifero, (drenaggi, evapotraspirazione, ricarica ecc) dipendenti in qualche modo dal carico idraulico rappresentano per ogni cella un termine dell’equazione (2). Per ogni cella si costruisce l’equazione (2) e quindi si risolve il sistema.
L’obiettivo di una simulazione transiente è generalmente quello di predire una distribuzione del carico idraulico in tempi successivi data una distribuzione iniziale di carico, i parametri idraulici e definiti gli stress esterni al sistema. La distribuzione iniziale del carico idraulico fornisce un valore h per ogni nodo della griglia all’inizio del primo intervallo di tempo (‘’time step’’). Il primo passo della soluzione è calcolare i valori di h al tempo t2 che marca la fine del primo time step per ogni nodo della griglia. Risolto il sistema di equazioni per il tempo t2 si passa al tempo t3 e così via per tutti gli intervalli di tempo fino a coprire il periodo di interesse. Ad ogni time step il sistema di equazioni viene riformulato.

2.4 Iterazione

Modflow utilizza un metodo numerico iterativo per ottenere la soluzione del sistema di equazioni implementate per ogni time step.
Il valore h da definire per la fine di un dato time step viene calcolato iniziando con l’assegnare un valore arbitrario o stima (ho), a ciascun nodo della griglia alla fine del time step stesso. Si comincia quindi una procedura di calcolo che varia queste stime, producendo un nuovo set di valori di h (h1)che sono più vicini alla soluzione del sistema di equazioni. Questi nuovi valori di h prendono il posto dei valori assunti inizialmente e si ripete il calcolo producendo un terzo set di valori (h2). Questa procedura si ripete producendo ogni volta un nuovo set di valori di h che soddisfano con approssimazione crescente il sistema di equazioni. Quando i valori di h approssimano quelli che soddisfano l’equazione le variazioni fra una iterazione e la successiva diventano molto piccole. Questo comportamento (convergenza) viene utilizzato per arrestare le iterazioni (criterio di chiusura). Esistono differenti metodi di iterazione: Modflow permette la scelta fra 4 di questi (vedi USGS Book 6-Modelling Technique). E’ importante notare che anche se ‘’esatte’’ le soluzioni di un modello sono sempre approssimazioni all’equazione differenziale di flusso (1).; l’accuratezza di questa approssimazione dipende da parecchi fattori, inclusi il criterio di chiusura scelto.


2.5 Tipi di celle e simulazione dei confini

Il sistema di equazioni sopra descritte non viene costruito per tutte le celle; Modflow raggruppa le celle escluse dalla computazione in due categorie ‘’celle a carico costante’’ e ‘’celle inattive’’ o celle ‘’no-flow’’. Le celle ad h costante sono quelle per cui il valore di h viene specificato all’inizio e resta costante attraverso tutta la simulazione; le celle inattive sono quelle per cui non è permesso alcun flusso in ingresso od in uscita. Le restanti celle della griglia vengono definite ‘’celle a carico variabile’’. Per ogni cella a carico variabile deve venir definito il sistema di equazioni di cui al paragrafo precedente. Celle inattive o a carico costante vengono in genere utilizzate per modellizzare i confini del sistema idrogeologico.


2.6 Discretizzazione verticale

Modflow discretizza lo spazio nel piano orizzontale definendo un numero di righe, un numero di colonne e l’ampiezza di ciascuna riga e colonna. La discretizzazione dello spazio verticale viene ottenuto specificando il numero di layers utilizzati. Modflow utilizza due differenti tecniche di discretizzazione verticale:

Ciascuno dei metodi di implementazione della discretizzazione verticale presenta vantaggi e difficoltà. Le equazioni del modello sono basate sull’assunzione che le proprietà idrauliche sono uniformi in ciascuna cella o che almeno parametri medi possono essere valutati per ogni singola cella; queste condizioni si incontrano quando il modello si conforma alle unità idrogeologiche (punto2). Inoltre ci si deve aspettare una maggiore accuratezza quando gli strati del modello corrispondono ad intervalli con perdite di carico trascurabili e anche questo è più probabile sotto la configurazione del punto 2. D’altra parte la griglia deformata non si adatta alle assunzioni su cui si basa il modello per esempio che celle individuali non dovrebbero avere facce rettangolari, che gli assi maggiori della conduttività idraulica siano allineati con gli assi del modello. Poiché ciascuno dei due metodi comporta degli errori, per minimizzarli la discretizzazione verticale viene ottenuta in pratica da unacombinazione dei punti 1. e 2.

3. APPLICAZIONE DEL MODELLO ALLA PIANA COSTIERA FRA IL TORRENTE TRIPESCE ED IL FOSSO DELLA MADONNA

MODFLOW è stato applicato alla piana costiera del Fiume Cecina da poco a monte della soglia idrogeologica della Steccaia fino al mare e dal T. Tripesce al Fosso della Madonna. L’applicazione del modello ad elementi finiti alla pianura costiera fra il torrente Tripesce ed il fosso della Madonna ha permesso la quantificazione del bilancio idrgeologico; quantificazione impossibile da realizzare con il metodo di bilancio classico in quanto risultavano sconosciuti parte degli afflussi nonché i deflussi a mare delle acque sotterranee. L’elevata densità di dati sia stratigrafici, idrogeologici e freatimetrici, raccolti da fonti bibliografiche e direttamente misurati (Fig.1) ha permesso l'elaborazione del modello idrogeologico concettuale e l’implementazione del modello numerico. In Fig. 3 sono riportate le ricostruzioni stratigrafiche effettuate sulla base dei pozzi lungo delle sezioni di Fig. 2 mentre in Fig 4 è riportata la stratigrafia idraulica come ottenuta da dati stratigrafici e prove di pompaggio. Da notare che in tutta la piana costiera risulta la presenza di due acquiferi principali spesso intercomunicanti localmente separati da pacchi di argilla di spessore talvolta rilevante. Le simulazioni calcolate hanno consentito la verifica della validità dell’interpretazione idrogeologica elaborata. Lo strumento informatico utilizzato, Visual Modflow, comprende infatti, un sistema di calibrazione per cui i dati ottenuti dalle simulazioni possono essere confrontati e valutati statisticamente con le misure freatimetrie reali effettuate su pozzi ‘’spia’’ o "calibration well" (non in pompaggio Fig.6- Fig.7); i dati di in ingresso del modello possono quindi essere corretti fino ad ottenere l’approssimazione richiesta della simulazione di taratura alla realtà fisicamente verificata.La delimitazione dell’area sopra descritta con celle inattive o a carico costante (mare) è risultata un compromesso fra le caratteristiche del modello (che si adatta con difficoltà a schematizzare realtà fisiche caratterizzate da rilevanti gradienti topografici), e l’interpretazione idrogeologica dell’area in studio che non può prescindere dal considerare Collemezzano e le Colline di Bibbona. In particolare la modellizzazione è stata eseguita in modo da ottenere due scenari finali simulati: il primo rappresentativo della fine della stagione secca, essenziale per calibrare, utilizzando i dati della campagna piezometrica autunnale, l’implementazione stessa ed il secondo per visualizzare la situazione a fine inverno.

3.1 Implementazione generale ed input

L’implementazione di un modello consiste nel:


3.1.1 Discretizzazione dello spazio e rapresentazione della superficie topografica
Il primo passo in una modellizzazione idrogeologica numerica è quello di discretizzare lo spazio fisico indagato nelle tre direzioni X,Y,Z.

La superficie topografica, immessa come terne di coordinate (Lat. Long, Quota s.l.m.; reticolo di riferimento Gauss-Boaga) si sovraimpone agli strati. Le quote, sia per la superfice che per la base degli strati vengono inserite come coordinate XYZ: MODFLOW provvede ad interpolare (metodo dell’inverso del quadrato della distanza) ricalcolando una quota per ogni cella. E’ chiaro che tanto più fitta è la discretizzazione sia verticale che orizzontale tanto più prossima alla realtà sarà la distribuzione delle quote interpolate.

3.1.2 Confini del modello e condizioni al contorno
In una modellizzazione di flusso di acque sotterranee è necessario definire e rappresentare i confini dell’area di studio nonché i corpi idrici superficiali quali mari, laghi, fiumi la iterazione con le acque sotterranee è tutt’altro che trascurabile.

3.1.2.1 Confini del modello
I confini del modello sono stati definiti sia utilizzando celle inattive che, poiché l’area di studio è un’area costiera, con celle a carico idraulico costante a rappresentare il mare. Le celle inattive sono state distribuite in modo da racchiudere l’area di modellizzazione: bassa pianura del Fiume Cecina dalla steccaia al mare e dal T. Tripesce al Fosso della Madonna. Lo strato 19 e parzialmente lo strato 18 sono stati considerati come limiti inferiori del modello e rappresentati con celle inattive.

3.1.2.2 Corpi idrici a carico idraulico costante
Il mare o un corpo idrico ad altezza presschè costante quale il laghetto della Steccaia sono modellizzabili in MODFLOW mediante apposita funzione GHB (General Head Boundary). Il flusso da una sorgente esterna, definita come GHB, in entrata od in uscita da una cella viene calcolato in proporzione alla differenza fra carico nella cella ed il carico assegnato alla sorgente esterna. I dati necessari per la definizione di una cella come GHB sono:

General head (GH): Altezza del livello dell’acqua sulla quota di riferimento
Conduttanza (Cg ): Parametro che rappresenta la resistenza al flusso fra un corpo idrico a livello costante (GHB) e le acque sotterranee

I valori di GH per la Steccaia sono stati ricavati dalla carta topografica 1:10.000; i valori di conduttanza sono stati calcolati in base a stratigrafie e prove di emungimento di pozzi dell’area.

3.1.2.3 Fiume
MODFLOW permette di considerare l’influenza di corpi idrici superficiali sulle acque sotterranee. I fiumi forniscono o drenano acqua dalle acque sotterranee dipendentemente dal gradiente idraulico esistente fra la superficie del fiume e le acque sotterranee stesse. Modflow simula i fiumi come corpi d’acqua separati dal sistema di acque sotterranee da strati di materiali a permeabilità costante, cella per cella. Per ogni cella del modello rappresentante il fiume sono richieste le seguenti informazioni:

Altezza del Fiume (H): quota s.l.m. del livello dell’acqua nel fiume per ogni periodo da modellizzare
Quota del fondo del Fiume (Ho): altezza sul livello del mare del letto del corso d’acqua.
Conduttanza (Cf ): parametro numerico che rappresenta la resistenza al flusso fra il corpo idrico superficiale ed l’acqua sotterranea.
La conduttanza si calcola con: lunghezza del ramo di fiume nella cella (L), ampiezza del fiume nella cella (W), spessore del letto dei sedimenti del fiume (M), conduttività idraulica del materiale del letto del fiume:

Cf = K L W / M

La modellizzazione del Fiume Cecina è stata uno dei passaggi più critici per l’implementazione del modello in quanto sono disponibili pochi dati sia per il calcolo della conduttanza sia per la definizione di H nel periodo estivo ed invernale (H1= Periodo estivo e H2 =Periodo invernale). Un tentativo per assegnare alle celle rappresentanti il fiume i parametri sopra indicati è stato quello di ricavare l’ altezza Ho dalla carta topografica 1:10.000 ed i valori di C ed H sia dalle informazioni di alcune prove di emungimento (ricavate da bibliografia) effettuate su pozzi in alveo, sia dalle informazioni ottenute per il calcolo del bilancio idrologico (idrometrogrammi). I dati così ottenuti sono stati inseriti nel modello e quindi corretti per quanto possibile in fase di calibrazione.

3.1.3 Proprietà delle celle
Ogni cella del modello non già definita come cella inattiva o cella a carico costante viene caratterizzata da un set di proprietà idrauliche:

L’esame di alcune prove di portata nonché l’insieme dei dati idraulici ricavate dalle domande di concessione dei pozzi ha permesso di ricavare le proprietà sopra indicate per le celle occupate dai pozzi. Queste proprietà sono state poi estrapolate sull’intera area in base all’interpretazione idrogeologica già elaborata. La calibrazione ha permesso di verificare ed affinare l’estrapolazione areale. Nella tabella che segue sono indicati i valori delle proprietà assegnate mentre nelle Fig 5c è rappresentata la distribuzione areale finale sullo strato del modello rappresentativo dell'acquifero più superficiale (Strato 10).

Proprietà

Kx = Ky (cm\sec)

Kz (cm\sec)

Ss (1\m)

Sy (-)

Por (-)

1

1.0 1.0 1.0e-09 0.99 0.99

2

1.0e-08 1.0e-08 0.1 1.0e-03 0.2

3

5.0e-05 5.0e-06 3.00e-03 0.2 0.25

4

1.0e-03 1.0e-03 1.0e-05 0.01 0.1

5

1.5e-03 2.5e-04 5.0e-05 0.05 0.08

6

0.05 0.05 1.0e-05 0.2 0.25

7 - 8

0.15 0.075 1.0e-06 0.3 0.35

9

5.0e-03 1.0e-03 2.0e-04 0.05 0.1

10

6.0e-03 1.5e-03 2.0e-04 5.0e-03 0.01

11

0.01 1.0e-03 5.0e-04 0.025 0.1

12

0.1 0.1 2.0e-04 0.2 0.3

La proprietà 1 attribuita alle celle che rappresentano il mare, celle a carico costante, è stata ‘’costruita’’ per simulare (compatibilmente con le limitazioni degli input imposte da VisualModflow) la capacità di ricezione e cessione idrica illimitate. Le proprietà da 2 a 12 sono state ricostruite dalle prove di emungimento, per classi di litotipi ragionevolmente omogenei sulla base delle stratigrafie disponibili (2-Argilla-acquicludo; 3-Livelli a matrice argilloso limosa-acquitardi; 4- sabbie limose; 5-Arenarie,tufi; 6-ghiaie sabbiose; 7-8- ghiaie sciolte; 9-sabbie mediamente compatte; 10 - intercalazioni di conglomerati ghiaie ed arenarie; 11-conglomerati; 12 - sabbie dunali).

3.1.4 Perturbazioni esterne al sistema
Tutto quell’insieme di fattori esterni al sistema idrogeologico ma in grado di condizionarne il comportamento e lo stato fisico reale vengono definiti come perturbazioni esterne. Tra essi si annoverano emungimenti, drenaggi, ricariche artificiali. Ricarica meteorica ed evapotraspirazione sono ‘’perturbazioni esterne’’ sempre presenti.



3.1.4.1 Pozzi

In MODFLOW la modellizzazione degli effetti dei pozzi (sia per emungimento che per ricarica) viene effettuata secondo le seguenti assunzioni:

dove, per un certo periodo di modellizzazione e per un particolare pozzo multistrato:Q1 è l’emungimento dallo strato 1; Qw è l’emungimento totale; T1 la trasmissività dello strato 1 e ST la somma delle trasmissività degli strati penetrati dal pozzo. Va sottolineato che una simile rappresentazione dei pozzi multistrato simula efficacemente gli effetti dell’emungimento. Per simulare l’effetto di by-pass idraulico costituito da pozzi multistrato abbandonati o inattivi è necessario ricorrere ad altre tecniche, comunque possibili con Modflow. Non avendo dati certi sulla distribuzione nel tempo degli emungimenti si è assunto in prima approssimazione che i pozzi ad uso idropotabile ed industriale fossero in massima parte attivi sia nel periodo estivo che invernale, mentre quelli ad uso irriguo fossero attivi solo nel periodo estivo. L’esame della carta delle differenze dei livelli piezometrici fra campagna primaverile ed autunnale (Fig 8) ha permesso di individuare quei pozzi per i quali l’ipotesi iniziale non risultava confermata e quindi variare l’implementazione. Per ciò che concerne l’entità degli emungimenti questi sono stati ricavati dalle dichiarazioni dei Proprietari registrate negli Uffici del Genio civile di Livorno.

3.1.4.2 Drenaggi
La funzione drenaggio è stata messa a punto per simulare gli effetti della rimozione di acqua da un acquifero con velocità proporzionale alla differenza di carico idraulico nell’acquifero e la quota dei dreni. La funzione si basa sull’assunzione che il drenaggio non ha effetto quando il livello dell’acqua sotterranea scende sotto la quota del dreno (h).

I parametri richiesti sono:
Quota di drenaggio (h): altezza del livello di acqua libera nel dreno.
Conduttanza (Cd): coefficiente che tiene conto delle perdite di carico fra il drenaggio e l’acquifero circostante.Nel caso della Val di Cecina la funzione drenaggio è stata utilizzata per simulare gli effetti della bonifica meccanica della Mazzanta.
Per l’intero periodo di simulazione si è considerato: h = -0.15 mslm e Cd = 500 m2/g
I valori sono stati ricavati da misure idrometriche effettuate nei canali di drenaggio.

3.1.4.3 Ricarica ed evapotraspirazione

 Ricarica

MODFLOW modellizza una ricarica areale dovuta alle precipitazioni, che percola attraverso gli strati fino alle acque sotterranee. I dati necessari per la modellizzazione sono stati ricavati dai calcoli del bilancio idrogeologico relativi alla piana costiera. In particolare si è considerata una infiltrazione efficace uniformemente distribuita su tutta l’area di studio pari a 53 mm/anno per il periodo estivo e pari a 180 mm/a per il periodo invernale.La ricarica viene applicata allo strato 1 ed il raggiungimento della tavola d’acqua viene calcolato da MODFLOW attraverso le proprietà idrauliche delle celle.

 Evapotraspirazione

La funzione Evapotraspirazione in MODFLOW simula l’effetto della traspirazione delle piante e la diretta evaporazione per rimozione di acqua dalle acque sotterranee. L’approccio è basato sulle seguenti assunzioni:

Come per i dati relativi alla ricarica anche i dati relativi all’evapotraspirazione sono stati ricavati dai dati di bilancio:
Ev Estate (giugno, luglio, agosto): 1720mm/a;
Ev Inverno (da settembre a maggio): 520 mm/a; Hex =1m

3.1.5 Pozzi per la calibrazione

VMODFLOW è dotato di un modulo specifico per consentire la calibrazione del modello implementato. La funzione calibrazione memorizza, per ogni periodo di modellizzazione, i valori di carico idraulico (Hc) calcolati per la cella ove sono ubicati i filtri di ciascun ‘’pozzo di osservazione’’. Ciò rende possibile paragonare i valori di carico idraulico calcolati con quelli osservati (Ho), valutare statisticamente la rispondenza del modello alla realtà, produrre idrogrammi per i pozzi di calibrazione. Per ogni pozzo di calibrazione vengono inseriti:

Coordinate X,Y del pozzo; Profondità del filtro; Periodo dell’osservazione; Altezza piezometrica osservata.

Per la modellizzazione in oggetto sono stati inseriti quei pozzi (circa 120) misurati durante le due campagne piezometriche di cui si conosceva almeno indicativamente la profondità del filtro. I dati della campagna piezometrica primaverile sono stati, come vedremo nel par. 3.1.7.1, utilizzati per definire la situazione iniziale; la calibrazione è stata effettuata utilizzando i dati piezometrici ricavati nella campagna autunnale.

3.1.6 Bilancio a zone

VMODFLOW permette il calcolo di un bilancio idrico per zone utilizzando i dati ottenuti dalla simulazione. Il dato di bilancio zona per zona viene fornito sia in totale che considerando separatamente gli elementi del sistema che contribuiscono agli scambi. Per ottenere questi dati di bilancio è necessario definire le zone di interesse. La Fig15 mostra le quattro zone distinte nella modellizzazione in oggetto.
Zona 1 Rappresentativa del mare . Tale zona è necessaria per controllare l’entità degli scambi con le zone adiacenti.
Zona 2 Area a Nord del Fiume Cecina
Zona 3 Area a Sud del Fiume Cecina
Zona 4 Porzione modellizzata a monte della Steccaia. Necessaria solo per completare il modello

3.1.7 Parametri per il calcolo

3.1.7.1 Regime transiente: situazione iniziale e discretizzazione del tempo

L’area di piana costiera modellizzata è come già ripetutamente accennato fortemente antropizzata e caratterizzata da forti emungimenti, variabili in intensità nel corso dell’anno, sia agricoli che industriali che idropotabili; ciò ha comportato l’implementazione e quindi il calcolo del modello per regime transiente. A differenza del calcolo per regime stazionario che prevede, data una situazione iniziale di carico piezometrico e considerate delle perturbazioni al sistema, il raggiungimento di uno stato di flusso stazionario il calcolo per regime transiente non prevedendo il raggiungimento dell’equilibrio permette la modellizzazione di un sistema idrogeologico sottoposto a stress esterni variabili nel tempo in numero ed intensità. Di tali stress esterni al sistema quali ricarica, emungimenti, evapotraspirazione si è già trattato nei paragrafi precedenti; è però necessario sottolineare che tali stress (soprattutto gli emungimenti) possono variare frequentemente in periodi brevi. La schematizzazione per cui essi vengono considerati costanti (per un periodo di 90 giorni, rappresentativo del periodo estivo, e di 275 giorni rappresentativo del periodo invernale) risulta una semplificazione indispensabile considerata l’incertezza sull’effettiva entità e sulla distribuzione temporale degli emungimenti. Il modello è stato quindi calcolato per due periodi ‘’ Stress period’’ uno di 90 giorni e l’altro di 275; per una maggiore accuratezza del calcolo matematico i due periodi sono stati a loro volta suddivisi in 10 intervalli ciascuno (‘’Time steps’’) all’interno dei quali si procede per l’iterazione come descritta al paragrafo 2.4. Il calcolo per regime transiente richiede, oltre alla discretizzazione del tempo l’implementazione di una situazione di carico idraulico di partenza. Come situazione iniziale sono stati utilizzati i dati relativi alle misure piezometriche della campagna primaverile.

3.1.7.2 Parametri generali
Oltre ai parametri sopra indicati che definiscono il sistema idrogeologico è necessario inserire alcuni dati necessari al calcolo.

  Caratteristiche idrologiche degli strati.
MODFLOW permette la definizione di quattro tipi di strato:
TIPO 0: Confinato -Trasmissività e coefficiente di immagazzinamento sono costanti per tutta la modellizzazione.
TIPO 1: Non confinato - La trasmissività dello strato varia. Viene calcolata dalla conduttività e dallo spessore saturo. Il coefficiente di immagazzinamento è costante.
TIPO 2: Confinato/Non confinato - La trasmissività dello strato è costante. Il coefficiente di immagazzinamento varia fra i valori per uno strato confinato e quelli per un non confinato.
TIPO 3: Confinato/Non confinato - La trasmissività dello strato varia. Viene calcolata dalla conduttività e dallo spessore saturo. Il coefficiente di immagazzinamento varia fra i valori per uno strato confinato e quelli per un non confinato. Gli scambi con gli acquiferi sottostanti sono limitati all’acquifero completamente desaturato.

Nel modello in esame lo strato 1 è stato definito come tipo 1, gli strati da 2 a 9 come tipo 2 e quelli da 10 a 18 come tipo 0.

 Caratteristiche per il bagnamento delle celle:
MODFLOW considera differenti modalità di bagnamento di celle secche; tale operazione richiede di definire:

  1. se il bagnamento avviene solo da celle sottostanti (piezometrica in risalita) o anche da celle laterali (caso più generale);
  2. il Fattore di bagnamento(W), questo parametro controlla il carico piezometrico (H) nella cella . Il carico piezometrico H è definito come:

H = BOT +W * (Hn-BOT)

dove: BOT = altezza della base della cella W = Fattore di bagnamento; Hn = altezza piezometrica delle celle circostanti che causa il bagnamento della cella in oggetto.
Per il caso in esame si è scelto di considerare possibile il bagnamento da tutte le celle circostanti e si è scello un fattore di bagnamento pari a 0.5.

 Iterazione: Il compito dell’iterazione è già stato descritto nel paragrafo 2.4. dove si è accennato all’esistenza di più metodi iterativi. Il metodo utilizzato in questo lavoro è denominato SIP ‘’Stronghly implicit procedure’’ ed è stato scelto perché nel caso specifico si è dimostrato più rapidamente convergente degli altri provati. Il massimo numero di iterazioni da eseguire è stato fissato a 100 per ogni time step. Tale numero è necessario solo per fermare il calcolo in caso di mancata convergenza.
Poiché ad ogni iterazione il programma controlla la massima variazione della soluzione (rispetto all’iterazione precedente) per ogni cella; se il massimo cambiamento è al disotto di una tolleranza preventivamente definita il programma passa ad un nuovo time step. Tale tolleranza o criterio di convergenza è stato imposto a 0.01 m.

3.2 Risultati della simulazione

Come già descritto nei paragrafi precedenti la simulazione è stata effettuata per ricostruire la distribuzione del flusso di acque sotterranee a fine estate (Time stress =1) ed a fine inverno (Time stress= 2) e per verificare la validità delle interpretazioni idrogeologiche implementate (per mezzo della calibrazione). Visual Modflow presenta la distribuzione dei vari parametri calcolati per ogni strato modellizzato.

Di seguito si illustrano i dati relativi a due strati rappresentativi degli acquiferi principali (Strato 10 e Strato 14) e si riportano le carte ed i dati relativi allo strato rappresentativo dell’acquifero superiore, lo strato 10. Poiché per la prima simulazione (Stress period =1) esistono dati piezometrici reali sia per l’inizio del periodo che per la fine, la funzione calibrazione può essere utilizzata nel pieno delle sue potenzialità


3.2.1 Carico idraulico
Confrontando le superfici piezometriche dei due strati considerati rappresentativi dei due acquiferi, riferite a ciascun periodo di simulazione, si nota che esse differiscono solo in minuti dettagli locali, mentre l’andamento generale appare pressochè identico; questa osservazione suggerisce che la separazione dei due acquiferi è efficace solo localmente, mentre a scala più vasta il sistema si comporta come un acquifero unico multistrato (La Fig 10 mostra il risultato della simulazione della piezometrica per lo strato 10 alla fine del periodo 1). Ciò appare confermato dalla elevata somiglianza della pezometrica "virtuale" (estiva) a quelle misurata, a dispetto della estrema variabilità delle profondità dei filtri dei punti di misura reale. La concordanza sostanziale fra piezometrica calcolata e quella osservata, come messo in evidenza dal grafico di calibrazione (Fig. 9), suggerisce una buona prossimità del modello interpretativo con la distribuzione reale delle caratteristiche idrauliche del sottosuolo e quindi l’attendibilità di ulteriori simulazioni. La simulazione relativa al Periodo 2 (inverno) è priva di termine di paragone; si sottolinea a tale proposito che mancando informazioni dettagliate circa l’effettiva distribuzione spaziale e temporale degli emungimenti per il periodo in questione; permane un margine di incertezza che può essere ridotto solamente acquisendo dati certi sugli emungimenti e/o eseguendo una nuova campagna di misure piezometriche a fine inverno che permetterebbe di identificare ubicazione ed entità effettiva degli emungimenti.


3.2.2 Velocità e direzioni di flusso
L’analisi dei risultati mostra che la direzione di flusso risulta pressochè identica per i due livelli acquiferi nei due periodi modellizati e che c’è una generale tendenza dell’alveo del fiume Cecina ad agire prevalentemente come drenaggio.

Alla luce di quanto esposto gli apporti del fiume alla falda sono da ascriversi soltanto come relativi ai periodi di morbida e piena.La carta della direzione di flusso, nello strato 10 considerato rappresentativo dell’acquifero superficiale, nell’area fra il Cecina e la Mazzanta mostra una evidente ingressione di acque marine verso l’entroterra. (La Fig 11 mostra la direzione del flusso nell'area della Mazzanta la Fig 12 La carta della conducibilità elettrica specifica per la stessa area). Per il periodo 1 in Fig 13 vengono riportate le velocità di flusso per lo strato 10. La lunghezza delle frecce è proporzionale al modulo del vettore velocità; nella tabella che segue sono indicate le velocità massime calcolate.Tali velocità massime si concentrano, in accordo con i marcati gradienti idraulici della zona, nell’area di Collemezzano; se ne deduce che, stante l’attuale distribuzione degli emungienti, anche ad uso idropotabile, prevalentemente al piede di Collemezzano, la prevenzione di eventuali sversamenti localizzati o inquinamenti diffusi è prioritaria quest’area.

3.2.3 Calibrazione

La Fig. 9 mostra il grafico di correlazione delle piezometrie calcolate con quelle osservate per tutti i pozzi di osservazione dell’area. Il grafico fornisce una indicazione del grado di approssimazione del modello alla realtà: in caso di corrispondenza esatta modello simulato- situazione reale tutti i punti dovrebbero giacere sulla retta a 45°. Per l’insieme dei pozzi di calibrazione sono stati calcolati l’errore medio, l’errore assoluto e l’errore quadratico medio.

Nella tabella che segue vengono mostrati i parametri sopra descritti calcolati oltre che per l’insieme dei pozzi di monitoraggio (Area tot) anche per sottogruppi:

Il risultato generale della simulazione risulta ampiamente soddisfacente; l’accuratezza del modello può essere definita:

  Area 0 Area 1 Area 2 Area 3 Area 4 Area 5 Area 6 Area Tot
Em

0.37

-0.21

1.05

-0.15

0.38

0.23

-0.30

0.09

Ea

0.80

0.60

2.13

0.23

0.73

0.76

1.26

1.02

RMS

0.97

0.66

2.55

0.31

0.84

0.89

1.59

1.44

Per migliorare l’accuratezza dei risultati è necessario in generale ottenere dati quantitativi più dettagliati sugli emungimenti nonché sulle caratteristiche dei pozzi di calibrazione; in particolare appare vivamente suggeribile l’implementazione di una modellizzazione a se stante della zona di Collemezzano, provvedendo ad un infittimento della griglia ed alla ricerca di dati il più possibile dettagliati sulla distribuzione temporale e sull’effettivo ammontare degli emungimenti, ed affinando la rappresentazione spaziale della distribuzione delle caratteristiche idrauliche del sottosuolo con le statigrafie di ulteriori perforazioni eventualmente eseguite in zona.

3.2.4 Idrogrammi
Le Fig. 14a e 14b mostrano il grafico Livelli piezometrici in funzione del Tempo per alcuni pozzi di osservazione situati per la Fig 14a nella zona di piana costiera fra la Mazzanta ed il Fiume Cecina e per la Fig 14b a Sud del Fiume Cecina. Il tempo corrisponde ai due periodi modellizzati 90 e 275 giorni rispettivamente (tempo totale 365 giorni).Mentre nella simulazione del periodo 1 si sono considerati in emungimento tutti i pozzi sia irrigui che industriali che idropotabili nel periodo sono stati disattivati tutti i pozzi irrigui e si sono ridotti gli emungimenti di quelli ad uso potabile o industriale. Nella Fig. 14a si può notare che a fine inverno i pozzi della piana costiera a Nord del F. Cecina non riacquistano completamente il carico idraulico iniziale cosa che invece avviene per i pozzi della piana costiera situati a Sud del Cecina; ciò conferma (ferme restando le limitazioni precedentemente esposte riguardo alla simulazione relativa al periodo invernale) quanto già dedotto con i metodi classici di bilancio idrogeologico: la piana costiera compresa tra le foci del Fiume Cecina e del Torrente Tripesce si trova in situazione di deficit di acque dolci.


3.2.5 Bilancio
I risultati finali ottenuti dal bilancio idrico sotterraneo per la piana costiera sono riportati in Fig.15. Per ciascuna zona e per ciascun periodo di simulazione sono stati calcolati:
Variazione dell’immagazzinamento;
Immagazzinamento;
Emungimenti da pozzi;
Drenaggi;
Ricarica meteorica;
Evapotraspirazione;
Scambi con il fiume;
Scambi con copi idrici a livello costante (GHB);
ottenendo quindi l’entità degli scambi sotterranei di ciascuna zona con le limitrofe. I dati più significativi si possono così sintetizzare:

Il modello implementato concorda nelle sue simulazioni con le osservazioni effettuate sul terreno. Il suo impiego per la verifica preventiva degli effetti di breve e lungo termine di futuri interventi sul territorio concernenti il prelievo o la immissione di acque nel sottosuolo può, già nella forma attale dell’implementazione, supportare adeguatamente scelte tecniche e decisioni al riguardo. È comunque doveroso sottolineare che una ampia serie di lacune conoscitive non consente, per ora, una simulazione di dettaglio, in particolare delle aree già attualmente oggetto di emungimenti concentrati e nell’immediata prossimità del corso del Cecina.

3.2.6 Limiti della modellizzazione
La simulazione sin qui descritta è stata preminentemente finalizzata alla verifica ed all’ eventuale aggiustamento dello schema idrogeologico dell’area di studio e del relativo bilancio idrologico; in tale processo sono emersi una serie di problemi e limitazioni che vengono esposte di seguito.E’ importante sottolineare che la maggioranza delle difficoltà incontrate sono da ascriversi a carenza di dati, mentre le caratteristiche e le potenzialità di simulazione dell’ "utensile Modflow" sono risultate pressochè sistematicamente sovrabbondanti, nonostante la relativa complessità dell’ ambiente modellizzato.

Il corso del Cecina
Come descritto nel par. 3.1.2.3, per modellizzare gli scambi idrici Fiume-falda, per ciascuna cella è necessario impostare, per ciascun periodo di simulazione, sia la quota dell’alveo che quella del pelo libero dell’acqua; il primo parametro è stato ricavato dalle carte topografiche con il miglior dettaglio disponibile (carta tecnica regionale scala 1:10000) mentre per il secondo si è dovuto, in assenza di dati statisticamente validi, ricorrere a correlazioni empiriche ed estemporanee tra gli idrometrogrammi di Ponte di Monterufoli ed osservazioni dirette del livello in vari tratti del fiume in vari regimi di portata: entrambe le scelte, anche se nel caso specifico inevitabili, riducono pesantemente l’accuratezza delle valutazioni dei rapporti fiume-falda. A tale proposito è opportuno ricordare che l’implementazione di un livello medio per un intero periodo di simulazione non consente di modellizzare adeguatamente gli effetti idrogeologici degli eventi di piena. Tale inconveniente è ovviabile solamente con una discretizzazione temporale molto fitta, con conseguente proporzionale espansione dei tempi di calcolo. La soluzione migliore resta in questo caso l’esecuzione di una serie di misure di portata istantanea, eseguite su sezioni significative dell’alveo in differenti regimi di portata, in base alle quali calibrare la simulazione per l’intorno del fiume.

Topografia accidentata
La rappresentazione della superficie topografica in un modello ad elementi finiti avviene mediante una ‘’gradonatura’’ delle pendici dei colli corrispondenti alle discretizzazioni in X,Y,Z. Le morfologie particolarmente accidentate nel caso specifico in particolar modo valli strette ed incise, vengono mediate dalla rappresentazione modellistica introducendo un errore sulle quote in corrispondenza e nell’immediato intorno delle valli stesse. Per ovviare a tale inconveniente, che porta ad evidenti discrepanze di calibrazione per i punti di misura lungo le valli, sarebbe necessario un infittimento della griglia X,Y che avrebbe però come conseguenza anche un marcato incremento dei tempi di calcolo. La griglia implementata nel caso specifico, come precedentemente accennato, rappresenta un compromesso fra accuratezza di calibrazione e scopi della modellizzazione (verifica a vasta scala). La possibilità di risolvere nel dettaglio l’idrogeologia locale di valli incise o altri marcati accidenti morfologici non viene comunque esclusa, essendo sempre possibile implementare una griglia fitta di dettaglio per le aree di interesse utilizzando per il modello a piccola scala, come condizioni al contorno l’output del modello generale.

Distribuzione temporale degli emungimenti
La discretizzazione temporale, per simulazioni di lungo periodo, è mirata prevalentemente ad identificare unità temporali omogenee per quanto riguarda ricarica, evapotraspirazione ed emungimenti. Nel caso specifico la discretizzazione temporale è stata definita basandosi principalmente sulla rilevante quantità di pozzi irrigui la cui attivazione si restringe in uno spazio temporale ben definito e prevedibile. Purtroppo le rilevanti utenze industriali ed idropotabili che insistono sull’area non hanno periodi di attivazione e spegnimento altrettanto regolari e preventivamente modellizzabili. Per simulazioni generali di area ciò non comporta che una piccola discrepanza fra dati misurati ed osservati, purchè le quantità totali di acque emunte, possano essere note. Ben diverso è il caso di simulazioni preventive di dettaglio mirate allo studio degli effetti e delle problematiche idrogeologiche locali relative all’attivazione di nuovi emungimenti cospicui; per implementare una simulazione accurata si necessita di una discretizzazione temporale molto fitta all’interno della quale sia possibile discriminare periodi di attivazione e spegnimento degli emungimenti più cospicui anche con cadenza inferiore al giorno.

4. CONCLUSIONI

In generale si può concludere che:

  1. Il modello implementato si è rivelato uno strumento di rappresentazione e calcolo sufficientemente versatile ed accurato. La verifica dello schema idrogeologico dell’area di studio, elaborato sulla base di dati diretti ed indiretti, è stata effettuata mediante la simulazione dell’evoluzione della situazione dalla fine del periodo invernale alla fine del periodo estivo estivo ed è stata confermata dalla buona correlazione fra dati osservati e calcolati.
  2. Il modello può, già nell’attuale implementazione, supportare scelte tecniche e decisioni consentendo la verifica preventiva degli effetti di breve e lungo termine di interventi sul territorio concernenti il prelievo o la immissione di acque nel sottosuolo.
  3. Mentre per la pianura costiera il modello ha già un alto grado di accuratezza, per ottenere simulazioni in grado di modellizzare nel dettaglio l’andamento della circolazione idrica sotterranea in zone a morfologia accidentata (colline di Bibbona) e sottoposte ad intenso emungimento (Collemezzano), è necessario (ove si preveda di intervenire) prevedere un’affinamento locale.

In dettaglio si osserva:

  1. la presenza di marcato stress delle falde in una zona di forma grossolanamente triangolare con vertice sul Fiume Cecina all’altezza di Podere S. Girolamo e base lungo la costa dalla foce del T. Tripesce ad 1 Km a sud della Foce del Fosso Vallescaia;
  2. la porzione dell’area di stress posta a nord del Fiume Cecina mostra una ingressione salina evidente fino a circa 2 Km dalla costa marcata per tutta l’area di bonifica della Mazzanta per entrambi gli acquiferi principali dell’area vista anche la ridotta efficacia della separazione idraulica fra i due.
  3. il flusso sotterraneo di alimentazione degli acquiferi della piana costiera proviene per la zona a nord del Cecina prevalentemente da Collemezzano con velocità di flusso piuttosto elevate;
  4. il flusso sotterraneo di alimentazione degli acquiferi della piana costiera per la zona a Sud proviene sia dalle Macchie dei Pianacci che dalle colline di Bibbona; non si rilevano zone di flusso concentrato e le velocità di flusso risultano moderate ed uniformi;
  5. è presente un flusso trasverso che, richiamato dagli ingenti emungimenti sottopassa il Cecina trasferendo le acque sotterranee della zona meridionale prospiciente il fiume stesso verso il piede di Collemezzano;
  6. l’alveo del Cecina dall’altezza delle Macchie dei Pianacci fino al mare funge per la maggioranza del tempo da dreno, contribuendo positivamente alla ricarica degli acquiferi solo durante le ondate di piene;
  7. la saltuaria ricarica del Fiume Cecina alle falde si sviluppa prevalentemente in destra idrografica, principalmente a causa della depressione della superficie freatica indotta in tale zona dai cospicui e diffusi emungimenti;
  8. dalla comparazione delle simulazioni estive ed invernali (in particolare dallo studio degli idrogrammi calcolati) si nota che la ricarica naturale ripristina i livelli della superficie freatica quasi ovunque nella zona a Sud del Cecina mentre ciò non avviene mai nella porzione fra il Cecina ed il Tripesce confermando quindi un evidente sovrasfruttamento;
  9. Il sovrasfruttamento sta conducendo ad una evidente riduzione della riserva per il concomitare di tre fattori in ordine di entità decrescente: emungimenti, effetto dreno dell’alveo del fiume e bonifica meccanica.

BIBLIOGRAFIA