DINAMICA LOCALE DEGLI INQUINANTI NELLE ACQUE DI FALDA: ADVEZIONE E

DISPERSIONE IDRODINAMICA

CAPITOLO 8

Acque sotterranee: la tipica stratificazione

- Zona insatura (o “zona vadosa”, vadose zone) - Zona satura di acqua - Falde acquifere (o “acquifero”, aquifer) - Rocce e materiali impenetrabili all’acqua (aquitard)

“My main point is that the water table should be defined only as thepressure surface where pore-water pressure is at local atmosphericpressure. Its definition should not refer to saturation. The top of thezone of saturation may be above, at, or below this surface. If adefinition for a nontechnical audience is required, the surfaceconnecting water levels in dug or shallow open wells will suffice inmost cases. In fact, it is instructive to note that the definition of thewater table as a pressure surface logically follows from how it ismeasured.”

Estratto tratto dall’articolo “The Water Table” di T. L. Holzer,Ground Water 48, 171 (2010)

Nota sul “Water Table”

Definizioni

volume di acqua:

volume totalew frazione (volumetrica) di acqua

volume di aria:

volume totalea frazione (volumetrica) di aria

Con riferimento ad un volume di campione di terreno:

w a n

w n nella zona satura

volume libero:

volume totalen porosità del materiale

pori/canali occupabili ad acqua/aria

Anche per la dinamica degli inquinanti trasportati dall’acqua nella matrice del terreno si adotta (come visto per la dinamica nell’aria e nelle acque superficiali) una equazione advezione-dispersione-reazione:

A ciò si aggiunge la diffusione molecolare (vera e propria) nel liquido.

La dispersione (rispetto al flusso medio) è principalmente di tipo meccanico, ed è dovuta alla tortuosità del percorso del liquido nella matrice del terreno (pori e canali di dimensioni diverse, ostacoli granulari di varie forme e dimensioni).

( , )h tD rMatrice di dispersione idrodinamica(il pedice ‘h’ sta per hydrodynamic)i cui elementi dipendono generalmentedal punto e dal tempo

Il processo complessivo di dispersione viene trattato come un processo“diffusivo” in un mezzo inomogeneo/anisotropo:

( , )( , ) ( , ) ( , ) ( , )h

c tt c t t c t

t

r

u r r D r r

Complessivamente, l’evoluzione temporale della concentrazione di inquinante in un punto r è descritta da

Eventuali fenomeni reattivi e/o di sink-source devono essere addizionati (per il momento li escludiamo).

Il primo passo è esplicitare il termine di advezione, cioè determinare l’evoluzione di ( , )tu r .

La zona satura (senza aria) e la zona insatura (con una frazione d’aria) hanno una diversa resistenza al trasporto di acqua per advection.

Procediamo trattando separatamente l’advezione nella zona satura e nella zona insatura.

Attenzione: la concentrazione ( , )c tr è qui espressa in termini di quantità di inquinante (es. massa) per unità di volume di terreno (cioè di spazio nel quale la dinamica ha luogo).

Advezione nella zona satura (100% acqua, 0% aria)

[Data la bassa velocità di scorrimento dell’acqua nel sottosuolo, vengonotrascurate le forze di attrito che si oppongono al moto.]

Nel sottosuolo, un elemento di acqua (centrato nel punto r e valutato altempo t) si muove a causa di forze agenti su di esso:

forza netta

La forza netta è data dalla risultante di varie forze:

forza peso dell’elemento stessogradiente del campo di energia gravitazionale

forza “di pressione”: è dovuta al fatto chela pressione è diversa in ogni punto dellasuperficie che delimita l’elemento

gradiente del campo di pressione

Vari contributi che concorrono a determinare la pressione in ognipunto del sottosuolo:

- peso della colonna di materiale sovrastante;

- depressioni per capillarità nei canali della matrice solida parzialmenteanidra (se siamo nella zona insatura);

- azioni esterne quali il pompaggio (estrazione) o l’immissione di acqua indeterminati punti: ciò fissa la pressione in tali punti.

zaumenta

aumenta

Il valore nullo è posto alla separazione tra zona satura e insatura.

I campi di energia potenziale e di pressione vengono “tradotti” in termini dicorrispondenti quote dette “carichi”:

( , )t r carico di pressione (pressure head):[Nella zona insatura è detto anche potenziale di matrice (matric potential)].Tale “quota efficace” dipende da tutte le coordinate del punto (non solo daz) e può variare nel tempo in quanto la pressione locale può cambiare.

carico di quota (elevation head) : zTale carico è associato al campo gravitazionale: è direttamente la quota zrispetto ad un arbitrario livello di riferimento (ad es. il livello del mare).

In condizioni stazionarie, cioè se ( , ) ( )t r r non varia nel tempo, per misurare ( ) r si può inserire un tubo (aperto sopra) perpendicolarmente al terreno fino al punto r; l’estremità inserita deve avere una membranapermeabile solo all’acqua; l’acqua salirà nel tubo fino ad una altezza pari a

( ) r .

La somma dei due contributi dà il carico idraulico (hydraulic head):

( , ) : ( , )h t z t r r

Posso determinare la velocità media del fluido ( , )tu r (che serve per il contributo di advection) dalla conoscenza del carico idraulico ( , )h tr ?

Direzione spontanea del flusso: verso punti a più basso carico idraulico

Legge di Darcy (1856)

/

/

/

x

y

z

operatoregradiente

conducibilità idraulica (hydraulic conductivity)[ La dimensione fisica di K è L T-1 (tipicamenteè espressa in cm/s o in m/giorno). ]

( , ) ( , )K

t h tn

u r r

Se il regime di moto dell’acqua è laminare (bassi numeri di Reynolds, indicativamente Re < 10), si verifica sperimentalmente che la velocità media è esprimibile secondo:

ˆ ˆ( , ) ( , )K

t h tn

d u r d r

La componente della velocità dell’acqua lungo una generica direzionespecificata da un versore si calcola medianted̂

proiezione sulla direzione scelta

ˆ ˆ( , ) ( , ) ( , ) ( , ) ( , )K K

u t t h t t h tn n

r u r r u r r

Se come direzione prendo proprio la direzione di scorrimento dell’acqua[cioè se considero ], tale operazione mi dà il modulo della velocitàdell’acqua:

ˆ ˆ ( , )td u r

Valore assoluto della derivata del carico idraulico nella direzione del flussomedio : tale quantità è denominata gradiente idraulico.

derivata del carico idraulico nella direzione del flusso medio (è negativa)

La relazione

( , ) ( , )K

t h tn

u r rstabilisce che:

1) il flusso istantaneo dell’acqua è diretto da punti ad elevato h a punti a più basso h (è costruita in modo tale che ciò sia soddisfatto);

(es. depressione creata da aspirazione in unpunto: “richiama” l’acqua circostante)

2) ( , )tu r è localmente perpendicolare alle superfici che congiungono

punti ad uguale valore di carico idraulico h (“superfici equipotenziali”).

La conducibilità idraulica K quantifica la resistenza del sistemacomplessivo (matrice + fluido) a “rispondere” dinamicamente algradiente idraulico.

l gK k

del mezzok permeabilità intrinsecaDimensione fisica L2 (tipicamente in cm2)

K dipende da:

- proprietà del mezzo poroso (dimensione e impaccamento delleparticelle);

- proprietà del liquido che fluisce (elevata densità incrementa iltrascinamento per gravità, bassa viscosità di shear riduce l’attrito efacilita lo scorrimento

La permeabilità intrinseca esprime la capacità di un dato mezzo poroso di lasciarsi attraversare da un qualsiasi fluido.

Empiricamente:2k C d

con C un “fattore di forma” (numero puro) dipendente dalla forma delle particelle del mezzo e dall’impaccamento, e d il diametro medio delle particelle.

In generale:

particelle di forma regolare,piccole dimensioni, alto impaccamento(es. argille)

bassa k

particelle di forma irregolare, grandi dimensioni, scarso impaccamento,(es. ghiaia)

elevata k

In termini di velocitàlineare, tipici valorisono 1 m/anno (si arrivaa 10 m/anno per terrenisabbiosi).

La dinamica è molto piùlenta rispetto a quella inacque superficiali e inaria!

Problema: in generale le condizioni non sono stazionarie!

È la situazione comune, ad esempio se si studia la dinamica delle acque sotterranee in seguito a pioggia intensa, o nel caso in cui si pompi acqua dal sottosuolo.

È necessario conoscere se/come ( , )h tr varia nel tempo: un cambio di carico idraulico può modificare in modo rilevante lo scorrimento dell’acqua sotterranea.

Dalla conoscenza di h(r,t) (suo valore in ogni punto al dato istante),mediante

posso determinare il campo di velocità attuale dell’acqua in ogni punto dellaregione di controllo!

( , ) ( , )K

t h tn

u r r

Occorre costruire un’equazione che descriva l’evoluzione di ( , )h tr

Consideriamo una celletta di volume V e valutiamo il tasso di variazione di volume d’acqua, /acquaV t , all’interno di essa:

Consideriamo i flussi di acqua entranti/uscenti attraverso le facce dellacelletta tenendo conto della porosità n del materiale:

con l’operatore divergenza.

( / 2, , , ) ( / 2, , , )

( , / 2, , ) ( , / 2, , )

( , , / 2, ) ( , , / 2, )

( / 2, , , ) ( / 2, , , )

( , / 2, , )

acquax x

y y

z z

x x

y y

Vu x x y z t u x x y z t n y z

tu x y y z t u x y y z t n x z

u x y z z t u x y z z t n x y

u x x y z t u x x y z tn x y z

xu x y y z t u

n

( , / 2, , )

( , , / 2, ) ( , , / 2, )

( , )( , ) ( , )( , )

z z

yx z

x y y z ty x z

yu x y z z t u x y z z t

n z x yz

u tu t u tn V V n t

x y z

rr ru r

rapporti incrementali

aree delle superfici effettivamente offerte al flusso di acqua

Sostituendo al posto di ( , )tu r la relazione si ha: ( , ) ( , )K

t h tn

u r r

2( , ) ( , )acquaV KV n h t V K h t

t n

r r

2 ( , )acquaVK h t

V t

r (§)

Dato il materiale, Ss rappresenta la variazione del volume di acqua immagazzinata al suo interno, per unità di volume di materiale, quando si applica una variazione unitaria di carico idraulico. E’ un parametro misurabile.

Per i terreni acquiferi nella regione satura si definisce il coefficiente di specific storage (immagazzinamento specifico) del materiale poroso:

: ( del mezzo)acquas

VS specific storage

V h

2 2 22

2 2 2x y z

operatore Laplaciano

Si ricava che

( , )acquas s

V h h tS S

V t t t

r (§§)

dal rapporto incrementalealla derivata

Uguagliando (§) e (§§) otteniamo l’equazione di evoluzione temporale cercata:

2( , )( , )

s

h t Kh t

t S

r

r

Il sistema di equazioni completo per ( , )h tr e ( , )tu r è il seguente:

2( , )( , )

( , ) ( , )

s

h t Kh t

t SK

t h tn

rr

u r r

equazione autonomaper ( , )h tr

La soluzione formale della prima equazione è semplice! [Basta notare che ha la stessa forma della equazione di diffusione in un mezzo omogeneo/isotropo, e ricordare il metodo di espansione in componenti di Fourier…].

La difficoltà maggiore consiste nell’includere le condizioni al contorno specifiche del caso che si sta trattando:

- specificare eventuali barriere confinanti e/o interfacce che consentono scambio di acqua con altri comparti; le condizioni di questo tipo possono essere anche molto complesse, ad esempio se sono presenti piante che assorbono acqua dal suolo e la disperdono nell’aria (sono termini di sink ?);

- specificare le condizioni iniziali: il profilo ( ,0)h r ;

- specificare l’eventuale controllo del carico idraulico, in determinati punti, da parte di azioni esterne. Ad esempio, azioni di pompaggio di acqua da punti specifici impongono una pressione locale specifica per garantire la portata di pompaggio stabilita; ciò impone che ( , )h tr abbia un valore stabilito in certi punti e sia controllato nel tempo.

Esempio di soluzione per un caso di interesse: pompaggio di acqua da un aquifer ideale confinato verticalmente, di spessore b, ma di estensione orizzontale infinita.

b

Il pozzo: aspira a portata costante in modouniforme lungo tutta l’estensione dellostrato confinato

pozzo di monitoraggio a distanza r

sezioni delle superficiisopotenziali (stesso carico idraulico): sono superfici cilindriche concentriche: il flusso medio di acqua è in orizzontale, radiale ed entrante nel pozzo…

stato iniziale: caricoidraulico costantein ogni punto

pompaggio 2 20( , ) 0.5772 ln( / 4 ) per / 4 0.1

4 s s

Qh r t h r S t r S t

K b

Per tempi t abbastanza lunghi dall’inizio del pompaggio e distanze piccole dal punto di pompaggio, una forma approssimata è

La soluzione analitica del problema è dovuta a C. V. Theis [Trans. Amer. Geophys. Union 16, 519 (1935)]:

pompaggio 20( , ) ( / 4 )

4 s

Qh r t h W r S t

K b

'

( ) ' per 0'

x

x

eW x dx x

x

pompaggioQ è la portata volumetrica dell’acqua pompata, 0h il carico idraulico

molto lontano dal punto di pompaggio, e ( )W x la “funzione pozzo” (è la “funzione esponenziale integrale”):

Vista dall’alto delle superfici equipotenziale (stesso carico idraulico):

h decrescente

Andando verso il centro la velocitàdell’acqua aumenta in quanto il gradienteidraulico cresce…

vettori a distanza fissadal centro (e ad un datoistante)

u

h

“Cono di drawdown”(ad un dato istante)

Spesso è importante riuscire a determinare in modo rapido (maapprossimativo) la velocità di scorrimento dell’acqua, in un dato punto,da poche misure dirette e senza fare calcoli complessi!

In condizioni stazionarie (carico idraulico ( )h r e velocità locale ( )u rche non variano nel tempo), per determinare la proiezione di ( )u r sul piano orizzontale (xy) è sufficiente misurare il carico idraulico in 3 punti nelle vicinanze del punto di interesse: è il “metodo dei tre punti”.

1) Misuro h1 , h2 , h3 in tre punti 1, 2 e 3 attorno al punto di interesse e posti alla stessa quota (tutti giacenti sullo stesso piano orizzontale); i punti siano etichettati in modo che 1 2 3h h h ;

3) Traccio la congiungente 1-3; assumo che il carico idraulico decresca linearmente andando da 1 a 3 e colloco un punto “A” su tale congiungente in corrispondenza del valore di h pari a h2;

2) Rappresento graficamente i tre punti (con le coordinate in scala effettiva);

4) Traccio la congiungente 2-A e la prolungo; considerando che 2h hA , assumo

che su tale retta tutti i punti abbiano lo stesso carico idraulico: è una approssimazione locale della intersezione della “superficie isopotenziale” (che contiene il punto 2) con il piano xy; la proiezione di u sul piano xy è quindi perpendicolare a tale retta e punta verso 3 (regione a basso carico idraulico);

5) Da 3 faccio partire la retta ortogonale alla congiungente 2-A; l’intersezione tra le due rette dà il punto B;

6) Valuto la proiezione di u sul piano xy dal gradiente idraulico tra il punto B (con 2Bh h noto) e il punto 3. Devo conoscere la conducibilità idraulica e la

porosità:

3 2 2 3orizz.

,3 ,3

nella regione tra 1, 2, 3B B

h h h hK Ku

n l n l

( ) ( )K h z z h z

n z

Per determinare la componente verticale della velocità media basta misurare il gradiente idraulico in tale direzione, cioè valutare il carico idraulico alla quota z e ad una quota z z :

1) lo scorrimento è nel verso della diminuzione di carico idraulico; 2) il modulo della componente verticale della velocità è pari a

Esempio. Per determinare la velocità media orizzontale dell’acqua in un punto della regione satura in condizioni stazionarie, si è misurato il carico idraulico in tre punti vicinali:

1

2

3

36.5 m35.0 m34.0 m

hhh

Noto che 610 m/sK e che la porosità del terreno è 0.32n , determinare direzione/verso/modulo della velocità media (sul piano orizzontale).

1,AlColloco il punto A sull’asse 1-3 tale che hA = h2 (assumo decrescita linearedel carico idraulico andando da 1 verso 3):

3 11 1,

1,3

11, 1,3

3 1

35 36.5 240 m 144 m

34 36.5

A A

AA

h hh h l

l

h hl l

h h

Traccio (graficamente) la retta che passaper 2 e per A; traccio la retta che passaper 3 e perpendicolare alla precedente;individuo B e misuro direttamente lalunghezza lB,3 . Poniamo che da questacostruzione risulti 85 m.

683 2

orizz.,3

10 m/s 35 34 3.7 10 m/s 3 mm/giorno

0.32 85B

h hKu

n l

Quindi: nei punti tra 1, 2, 3 l’acqua scorre nella direzione/verso B3,con velocità (in orizzontale) di 3 mm/giorno.

Con questa semplice tecnica possostabilire in che direzione e con chevelocità si sta propagando un inquinanteper advezione!

Advezione nella zona insatura (mix aria/acqua)

La dinamica di advezione nella zona insatura è analoga a quella vista per la zona satura. La differenza consiste nel fatto che sia la conducibilitàidraulica che il carico idraulico dipendono dal grado di idratazione della matrice, che è espresso dalla frazione volumetrica di acqua:

volume di acqua:

volume di terrenow

Ulteriore complicazione: in ogni punto, w può cambiare nel tempo!

La frazione di acqua nella zona insatura varia tra un valore minimo rdetto “contenuto di acqua irriducibile” e il valore massimo (sul water table) che è pari alla porosità della matrice: w n .

r (ad una data temperatura) corrisponde al contenuto minimo di acqua che non può essere estratto dalla matrice neppure creando elevatissime depressioni (costituisce il limite all’essiccamento del terreno).

Tipici profili sono i seguenti. Notare che le curve in “andata” (wetting) e in“ritorno” (drying) non coincidono: fenomeno di isteresi.

La conducibilità idraulica aumenta all’aumentare della frazione di acqua. Più la matrice è impregnata di acqua, più risulta permeabile all’acqua stessa. Si spiega con il fatto che i pori sono connessi da canali che contengono già dell’acqua; ciò facilita lo scorrimento come processo collettivo…

Il carico di pressione è negativo nella zona insatura (si annulla alla saturazione)

( )w

( )wK

[ ( , ) ( , )]acquaw

VV t t

t

r u r

Come già visto per la zona satura, valutiamo la variazione del volume diacqua (Vacqua) in una celletta di volume V, nel tempo t. Ripetendo losviluppo si arriva a

( / ) ( , )acqua wV V t

t t

rNotiamo che

( , )[ ( , )] ( , )w

w

tK t h t

t

r

r rSegue quindi

( , ) ( , )[ ( , )]

( , )[ ( , )]

( , )[ ( , )] 1

ww

w

w

z t h z tK z t

t z zz z t

K z tz z

z tK z t

z z

L’assunzione comunemente fatta per la zona insatura è di uniformità sulpiano xy; si considera solo la componente z:

( , ) [ ( , )]wz t z t

Dal fatto che si ha corrispondenza biunivoca tra carico di pressione efrazione d’acqua (per il wetting), cioè esiste una funzione , si puòesprimere

( )w

Per risolverla occorre utilizzare i profili empirici di wetting del tipomostrato sopra:

( )wK ( )w

( , )

( , ) ( , )[ ( , )] 1

w

w ww

w z t

z t z tdK z t

t z d z

Si arriva ad una equazione di evoluzione autonoma per :

E’ una forma particolare della equazione di Richardsper la zona insatura

( , )w z t

Occorre conoscere le condizioni iniziali, cioè il grado di idratazione in ogni punto della zona insatura al tempo-zero: ( ,0)w r

Ottenuta la soluzione ( , )w t r si procede in sequenza:

( , )w t r

( , )

( , )

( , ) ( )

( , ) ( )w

w

w t

w t

h t h

K t K

r

r

r

r

( , )( , ) ( , )

( , )w

K tt h t

t

ru r r

r

ho costruito il termine di advezione!

nella zona insatura entra lafrazione di acqua al postodella porosità n

Contributo di dispersione idrodinamicanelle zone satura/insatura

Ricordiamo che il contributo di dispersione dell’inquinante è modellizzato come un processo di “diffusione”:

dispersione

( , )( , ) ( , )h

c tt c t

t

r

D r r

Con riferimento ad uno specifico sistema di assi principali*, la matrice didispersione risulta diagonale:

,

,

,

( , )

( , ) ( , )

( , )

h x

h h y

h z

D t

t D t

D t

r

D r r

r

coefficienti di dispersione idrodinamicanelle direzioni principali

asse z (verticale): ortogonale al suoloasse y (longitudinale): in orizzontale e orientato nella direzione media di scorrimento dell’acquaasse x (trasversale): il rimanente asse

* Mancano indicazioni chiare su come individuare tali assi principali. Scelta ragionevole:

Empiricamente si riscontra che:

,

,

,

( , ) ( , )( , ) ( , )( , ) ( , )

h x x m

h y y m

h z z m

D t u t DD t u t DD t u t D

r rr rr r

I tre coefficienti di dispersione comprendono due contributi:

- la dispersione meccanica dovuta alle inomogeneità dei granuli della matrice (contributo proporzionale alla velocità media dell’acqua);

- diffusione molecolare nell’acqua all’interno dei pori (unico coeff. di diffusione mD )

Tipici valori del coefficienti di dispersione in acque sotterranee: tra 10-2 e10 cm2/s.

dispersività del mezzo nelle tredirezioni (dimensione fisica: L).[Si possono determinare sperimentalmentemonitorando la dispersione di “puffs” ditraccianti rilasciati nel terreno.]

Torniamo all’equazione advezione-dispersione…

( , )( , ) ( , ) ( , ) ( , )h

c tt c t t c t

t

r

u r r D r r

Abbiamo discusso come ottenere la velocità media dell’acqua, ( , )tu r , nelle zone satura/insatura, e come quantificare i coefficienti di dispersione.

Estensioni da introdurre (se rilevanti, da valutare caso per caso):

1) In presenza di immissione diretta (source) di inquinante dal comparto-terreno occorre aggiungere uno specifico termine ( , )S tr al 2° membro.

- Se il rilascio di inquinante è istantaneo al tempo-zero, basta solo definire la specifica condizione iniziale ( ,0)c r che viene creata (non c’è un termine

( , )S tr da aggiungere…).

- Se invece il rilascio si protrae nel tempo (ad esempio se l’inquinante fuoriesce da una tubatura rotta…) allora occorre inserire il corrispondente termine ( , )S tr .

immissione da t = 0

profilo a t = 0

no ( , )S tr occorre ( , )S tr

( ,0)c r

2) In presenza di reazioni chimiche occorre aggiungere l’appropriato termine [ ( , )]R c tr .

3) Trattare l’eventuale ripartizione dell’inquinante tra acqua e materiale della matrice solida (per assorbimento/adsorbimento).

4) Se l’inquinante è volatile, nella regione insatura occorre trattare anche la ripartizione dell’inquinante tra acqua ed aria.

Nel seguito sviluppiamo i punti 2), 3), 4).

(Alcune) ulteriori complicazioni:

a) Anche l’aria, nella regione insatura, ha una propria dinamica e può trasporare l’inquinante anche all’esterno del comparto-terreno.

b) Le cinetiche di scambio di inquinante acqua-solido e acqua-aria possono essere lente! Ciò complica l’analisi…

c) Nel comparto-terreno, oltre all’acqua, possono essere presenti anche altri liquidi non solubili in essa, ma nei quali l’inquinante è solubile. Anche tali liquidi sono soggetti ad una propria dinamica!

Effetto “ritardante” della ripartizione dell’inquinante, tra acqua e matrice solida, sulla dinamica (regione satura)

Introduciamo il coefficiente di ripartizione della specie soil-water in condizioni di equilibrio a temperatura fissata:

,,

,

S eqS W

W eq

CK

c

concentrazione volumetrica di inquinante in fase acquosa: massa inquinate / volume di acqua

concentrazione massa/massa di inquinante nella matrice solida: massa inquinante / massa di solido

Assumiamo che i processi cinetici di scambio di inquinante tra fase acquosa e solido siano estremamente rapidi: equilibrio di ripartizione ad ogni istante e in ogni punto della regione. Quindi

,

( , )

( , )S

S WW

C tK

c t

r

r ,( , ) ( , )S S W WC t K c tr r

(come nel “Livello II” di Mackay …)

Consideriamo una celletta di materiale, di volume V . Poniamo:

totM la massa totale di inquinante presente nella celletta (parte in acqua e parte nel solido) al tempo t ;

SM la massa di solido presente nella celletta;

S “bulk density” del terreno = massa di solido / volume totale (acqua inclusa).

tot , in al tempo ( , )S W S WM V t V K n c t r

Bilancio di massa rispetto all’inquinante:

tot

,

,

in al tempo ( , )( , ) ( , )( , ) ( , )

S S W

S W W S W

S W W S W

M V t C M c t n VK c t M c t n VK c t V c t n V

rr rr r

volume di acquanella celletta

,tot,

( , ) ( , ) in S W SW

S W S

V K nM c t c tV V K n

t t n t

r r

Passiamo al tasso di variazione di massa in V:( , ) ( , ) /Wc t c t nr r

D’altro canto, la massa di inquinante dentro la celletta può variare solo peradvezione e dispersione (per ora ignoriamo reazioni e source/sink).

Uguagliando le due espressioni si ottiene

, ( , )( , ) ( , )

( , ) ( , )

S W S

h

K n c tV V t c t

n tV t c t

ru r r

D r r

Dividendo membro a membro per V otteniamo l’equazione finale:

( , )( , ) ( , ) ( , ) ( , )h

c tt c t t c t

t

r

u r r D r rR

,1 S W SK

n

R “fattore di ritardo”

(retardation factor)

tot in ( , ) ( , )

( , ) ( , )h

MV V t c t

tV t c t

u r r

D r r

Osserviamo che 1R ed è un numero puro tanto più elevato

- quanto più la densità del terreno è alta; - quanto più l’inquinante è “affine” solido; - quanto più la porosità è bassa.

E’ l’usuale eq. di advezione-dispersione, ma è come se l’acqua avesse la velocità fittizia '( , )tu r con modulo inferiore a quella vera ( '( , ) ( , )u t u tr rin quanto 1R ).

Per questo motivo si usa il termine “ritardo”: di fatto, lo scambio di inquinante tra acqua e solido rallenta la sua dinamica di spostamento nello spazio senza però modificare il tipo di evoluzione [è una sorta di “effetto moviola”…].

( , )'( , ) ( , ) ' ( , ) ( , )h

c tt c t t c t

t

r

u r r D r r

( , )'( , )

tt

u ru r

R,

,

,

' ( , ) '( , )' ( , ) '( , )' ( , ) '( , )

h x x

h y y

h z z

D t u tD t u tD t u t

r rr rr r

(i contributi Dm di diffusionemolecolare sono stati trascurati)

con

Dividiamo membro a membro per R:

Più R è elevato, più l’inquinantesi sposta /disperde lentamente.

Dato il tipo di terreno (densità eporosità fissate), R cresce alcrescere della costante diripartizione Ks,w per la speciechimica in esame.

Inquinanti diversi si spostano/disperdono nel terreno convelocità diverse [osservarel’analogia con la cromatografiain colonna solido /liquido!]

Effetto “ritardante” della ripartizione dell’inquinante tra acqua, matrice solida, e aria (regione insatura)

La trattazione è identica a quella vista per la regione satura. Occorre introdurre anche la costante di ripartizione aria-acqua:

,,

,

A eqA W

W eq

cK

c

con ,W eqc la concentrazione volumetrica (massa/volume di acqua) di

inquinante nella fase acquosa e ,A eqc la concentrazione volumetrica

(massa/volume di aria) di inquinante nell’aria.

Dai bilanci di materia (farlo come utile esercizio!) si ottiene ancora una equazione advection-dispersione modificata :

( , )( , ) ( , ) ( , ) ( , )h

c tt c t t c t

t

r

u r r D r rR

in cui ora il fattore di ritardo è dato da:

, ,1 S W S A W a

w

K K

R

con a wn la frazione di aria nella regione insatura. [Per i passaggi tenere presente che nella regione insatura vale ( , ) ( , ) /W wc t c t r r ].

Assumiamo:

1) che gli scambi di inquinante tra acqua-aria-solido siano molto rapidi, tali in ogni istante si abbia equilibrio di ripartizione;

2) che nella regione insatura w sia la stessa in ogni punto e non vari nel tempo.

Effetto di reazioni di decomposizione del 1° ordine in acqua, aria, solido

Ammettiamo la possibilità che l’inquinante si decomponga in acqua, in aria e sul solido secondo le seguenti reazioni

(reazione in acqua) Inquinante prodotti

(reazione in aria) Inquinante prodotti

(reazione nel solido) Inquinante prodotti

w

a

s

k

k

k

con leggi cinetiche del 1° ordine così scritte:

reaz.

reaz.

reaz.

( , )( , )

( , )( , )

( , )( , )

Ww W

Aa A

Ss S

c tk c t

tc t

k c tt

C tk C t

t

rr

rr

rr

in acqua

in aria

nel solido

costanticinetiche

Dai bilanci di materia (farlo come utile esercizio!) si ottiene ora una equazione advezione-dispersione-reazione:

( , )( , ) ( , ) ( , ) ( , ) ( , )h eff

c tt c t t c t k c t

t

ru r r D r r r

R

con il fattore di ritardo (già visto)

, ,1 S W s A W a

w

K K

R

, ,s S W S a A W aeff w

w

k K k Kk k

con keff una “costante cinetica efficace“ data da

Esercizio. 1000 gr di un liquido organico fuoriescono da un pozzo che pesca nella zona satura del sottosuolo. Noto che: - 0.1 m/giornoK (conducibilità idraulica), - 0.33n (porosità), - il gradiente idraulico vale in modulo 1.1 nella direzione verticale e 0.1 nella direzione orizzontale del flusso, - 1 mx y z (coefficienti di dispersività),

- 5 21.02 10 cm /smD (coeff. di diffusione molecolare nell’acqua)

Determinare:

1) a quali coordinate si ha il picco della concentrazione volumetrica di inquinante dopo 30 giorni dalla fuoriuscita (prendere come coordinate0,0,0 quelle del punto di rilascio);

2) quanto vale la concentrazione a tali coordinate se l’inquinante si ripartisce tra acqua e matrice solida, noto che la densità “di bulk” del terreno è pari a 32 gr/cmS e il coefficiente di ripartizione vale

, 0.5S WK litri acqua per kg di terreno.

Caso 1)

Scegliamo il sistema di assi:

Valutiamo le componenti della velocità media:

-5

-4

00.1 m/giorno

0.1 0.030 m/giorno = 3.5 10 cm/s0.33

0.1 m/giorno1.1 0.333 m/giorno = 3.9 10 cm/s

0.33

x

y

z

uK dh

un dyK dh

un dz

Ricordiamo le relazioni per i coefficienti di dispersione idrodinamica:

,

,

,

h x x m

h y y m

h z z m

D u DD u DD u D

Il profilo di concentrazione di inquinante è una Gaussiana (in tre dimensioni) con allargamenti omogenei nelle tre direzioni Cartesiane, e il cui centro si sposta con la velocità media dell’acqua:

2 220

3/ 2, , ,, , ,

1( , , , ) exp

44

y z

h x h y h zh x h y h z

y u t z u tm xc x y z t

t D D Dt D D D

2 2 6 2, , , 3.92 10 cm /s = 3.92 10 m /sh x h y h zD D D

Dato che le tre dispersività sono uguali si ottiene che

1/ 22 2 2 43.92 10 cm/sx y zu u u u

Valutiamo gli spostamenti del picco di concentrazione dopo

30 giorni = 2.6 106 s

Lo spazio percorso è nullo nella direzione x, mentre nella direzione orizzontale del flusso e in quella verticale si ha:

-5 63.5 10 cm/s 2.6 10 s 91 cm (0.91 metri)yu t 4 63.9 10 cm/s 2.6 10 s 1014 cm (10.1 metri)zu t

La concentrazione massima (picco) vale

3 / 2 36 6

4 33

(0, 0.91 m, 10.1 m 30 giorni)1000 gr

4 2.6 10 s 3.92 10

1000 gr6.9 10 gr/dm di terreno

1450 m

c y z t

(attenzione: per ottenere la concentrazione di inquinante nell’acquapresente nel terreno devo dividere il valore ottenuto per la porosità n).

Caso 2)

3 3 3 3, 2 10 kg/m 0.5 10 m / kg

1 1 4.00.33

S S WK

n

R

La soluzione della equazione advection-diffusione è la stessa di prima, salvo il fatto che occorre usare velocità e coeff. di dispersione “corretti”per il fattore di ritardo:

2 2' '2

0' ' '3 / 2 ' ' ', , ,, , ,

1( , , , ) exp

44

y z

h x h y h zh x h y h z

y u t z u tm xc x y z t

t D D Dt D D D

-5 -6

-4 -5

' 0' / = 3.5 10 cm/s /4.0 = 8.8 10 cm/s

' / = 3.9 10 cm/s /4.0 = 9.8 10 cm/s

x

y y

z z

uu u

u u

RR

, , ,

' ' ', , ,= / = / = /

h x h y h zh x h y h zD D D D D DR , R , R

, , ,

' ' ' 6 2 7 2= 3.92 10 / 4.0 m /s = 9.80 10 m /sh x h y h z

D D D

Le nuove coordinate del picco di concentrazione sono:

-6 6' 8.8 10 cm/s 2.6 10 s 22.9 cm (0.23 metri)yu t -5 6' 9.8 10 cm/s 2.6 10 s 255 cm (2.6 metri)zu t

La concentrazione, nello stesso punto, è molto più bassa! L’inquinante ha “corso” meno a parità di tempo di attesa….

3 / 2 36 7

2 2

6 7 7

5.56 5 33

(0, 0.91 m, 10.1 m 30 giorni)1000 gr

4 2.6 10 s 9.80 10

0.91 0.23 10.1 2.61exp

4 2.6 10 9.80 10 9.80 10

1000 gr2.1 10 gr/dm di terreno

181 m

c y z t

e

Se valuto la concentrazione in corrispondenza della posizione del picco delcaso precedente ottengo