Teoria del modello LAPMOD

LAPMOD (Bellasio and Bianconi, 2012a; Bellasio et al. 2017; Bellasio et al. 2018, Bellasio and Bianconi, 2022; Bellasio and Bianconi, 2025; Bellasio and Bianconi, 2026) è un modello Lagrangiano a particelle, tridimensionale e non stazionario, multisorgente e multisostanza, adatto a simulare la dispersione in atmosfera di sostanze inerti o radioattive, emesse sia in fase gassosa sia in forma di aerosol. Oltre alla dispersione degli inquinanti convenzionali, LAPMOD è in grado di simulare anche la dispersione in atmosfera di sostanze odorigene. Le particelle utilizzate per il calcolo sono spostate nell’atmosfera per effetto del vento (avvezione) e della turbolenza atmosferica (dispersione). Ciascuna particella trasporta una frazione della massa di inquinante emessa. Ad ogni instante è possibile calcolare la concentrazione e la deposizione in corrispondenza dei recettori considerando la posizione delle particelle e le loro masse.

LAPMOD è completamente accoppiato con il modello meteorologico diagnostico CALMET, fino alla versione 5.8.5, il quale è in grado di fornire tutte le informazioni necessarie riguardo la velocità e la direzione del vento e i parametri turbolenti.

Emissioni

LAPMOD permette di simulare la dispersione in atmosfera di più sostanze contemporanemante, emesse da una o più sorgenti di diverso tipo.

Le sostanze emesse possono essere inerti o radioattive, sia in fase gassosa sia in forma di aerosol. Oltre alla dispersione degli inquinanti convenzionali, LAPMOD è in grado di simulare anche la dispersione in atmosfera di sostanze odorigene.

Il rateo di emissione è sempre espresso in termini di massa / tempo (o attività / tempo o unità odorimetriche / tempo), indipendentemente dalla tipologia di sorgente.

Tipologie di sorgente

• Sorgenti puntuali con o senza spinta termica
• Sorgenti lineari
• Sorgenti areali
• Sorgenti a box tridimensionale
• Sorgenti sferiche
• Sorgenti poligonali di forma arbitraria

Rilasci di sostanze radioattive

Il rateo di emissione, in caso di rilasci di sostanze radioattive, è espresso in Bq/s (o multiplo di Bq). Il modello LAPMOD trasforma il rateo di emissione di attività radioattiva in rateo di emissione di massa (cioé massa emessa nell'unità di tempo), che può essere calcolato attraverso l'attività specifica (Knoll,1979), essendo la massa pari al rapporto tra attività e attività specifica.

Indicando con `A` l'attività (in Bq/s), con `M` il peso molecolare, con `A_0` il numero di Avogadro, con `lambda` il coefficiente di decadimento (in s) vale la relazione:

`massa = (A M) / (lambda A_0)`

Tra la costante di decadimento `lamda` ed il tempo di dimezzamento `T_(1/2)` vale la relazione:

`T_(1/2) = log 2/ lambda`

Il fattore di conversione tra il rateo di emissione di attività radioattiva (in Bq/s) e il rateo di emissione di massa (in g/s) è dato dal fattore `M/(lambda A_0)`, che è una proprietà del radionuclide emesso. Naturalmente LAPMOD può essere usato anche per simulare la dispersione in atmosfera di sostanze inerti (vale a dire non radioattive). Per queste sostanze il fattore di conversione è semplicemente pari a 1.

La trasformazione del rateo di emissione da Bq/s a g/s viene effettuata all'interno della routine SETSUB.

Fase gassosa e fase aerosol

Le sostanze rilasciate in atmosfera possono essere sia in fase gassosa sia in fase aerosol.

Nel caso degli aerosol si possono simulare particelle di taglia singola e particelle con una distribuzione log-normale della massa, avente un AMAD (Aerodynamic Mean Aerosol Diameter, Diametro Aerodinamico Medio dell'Aerosol) e GSD (Geometrical Standard Deviation, Deviazione Standard Geometrica) che dipendono dalla sostanza e dal tipo di rilascio. La distribuzione log-normale continua viene rappresentata nel modello in modo discreto, suddividendola in un certo numero di intervalli (NBINS, in genere uguale a 10) ciascuno con la propria probabilità. Per ogni particella di aerosol emessa si generano quindi NBINS particelle computazionali, ciascuna con diametro individuato da uno degli NBINS intervalli in cui è stato suddiviso l'AMAD.

Anche la distribuzione in numero di particelle di aerosol è una distribuzione log-normale. La probabilità di trovare una particella con diametro inferiore a `l` (che matematicamente è rappresentata dall'integrale della curva della distribuzione con `x < l`) è data da:

`Prob( x < l ) = 0.5 (1+erf(a))`

Dove:

`a = (log(l) - log(m)) / (sqrt 2 log (sigma_g))`

e `m` è la media geometrica (che è la mediana di una distribuzione log-normale) e `sigma_g` è la deviazione standard geometrica.

La distribuzione in volume delle particelle è proporzionale alla distribuzione di massa delle particelle attraverso la densità (costante) ed è anch'essa una distribuzione log-normale (Seinfeld e Pandis, 1998, eq 7.51), con la stessa deviazione standard geometrica della distribuzione in numero.

Il diametro mediano per la distribuzione in volume è dato da (Seinfeld e Pandis, 1998, eq 7.52):

`log( mVol ) = log(m) +3log^2(sigma_g)`

Quindi, per la distribuzione in volume, la probabilitè di avere un volume inferiore a V risulta:

`Prob( nu < V ) = 0.5 (1+erf(A))`

dove:

`A = (log(l) - log(mvol)) / (sqrt 2 log (sigma_g)) = (log(l)-log(m)-3log^2(sigma_g)) / (sqrt 2 log (sigma_g))`

Il problema si riduce perciò a trovare il valore di `A` tale che `Prob( nu < V )` sia un multiplo di (1/NBINS), cioè trovare un valore di `A` che soddisfi:

`erf(A) = 2 Prob( nu < V ) -1 `

Noto `A` il corrispondente valore di `l` può essere calcolato attraverso la relazione:

`l = exp (A sqrt 2 log( sigma_g ) + log(m) +3log^2(sigma_g))`

Valori tipici di AMAD e GSD per alcuni radionuclidi sono proposti da Baklanov e Sorensen (2001), e riportati nella seguente tabella.

Nuclide	Tempo di dimezzamento (ore)	AMAD (`mu m`)	Velocità di deposizione (cm/s)	GSD (`mu m`)
Cs-137	2.628E05	0.68	0.1	1.8 - 2.5
Cs-134	1.805E04	0.59	0.12	1.8 - 2.5
I-131	1.930E02	gas / 0.48	0.6	3.0 - 4.0
Te-132	7.820E01	0.81	0.3	1.5 - 2.5
I-133	2.080E01	gas / 0.60	0.7	non disponibile
Ba-140	3.048E02	0.45	0.9	non disponibile
Sr-90	2.549E05	2.50	2	2.0 - 2.5
Ru-103	9.459E02	0.65	0.5	1.4 - 2.5
Pu-238	7.683E05	4.30	2	2.0
Pu-239	2.111E08	4.30	2	2.1
Sr-89	1.212E03	2.50	2	non disponibile

Si può notare che per l'isotopo I-131, durante la prima settimana dopo il rilascio, domina la fase gassosa.

Input meteorologico

L'input meteorologico di LAPMOD è costituito da campi tridimensionali di vento e temperatura e campi bidimensionali di parametri turbolenti quali lunghezza di Monin Obukhov, velocità di frizione, velocità di scala convettiva, altezza dello strato di rimescolamento.

LAPMOD può leggere direttamente i campi meteorologici generati in output da CALMET (Scire et al., 1999a).

CALMET è il modello meteorologico diagnostico usato anche per preparare l'input del modello di dispersione CALPUFF (Scire et al., 1999b), che appartiene alla lista di modelli suggeriti da US-EPA. Sia CALMET sia CALPUFF sono disponibili sul sito web di Exponent (http://www.src.com).

LAPMOD è completamente accoppiato con CALMET fino alla versione 5.8.5, la più recente approvata dell'EPA.

Anche le variabili geofisiche necessarie a LAPMOD derivate dalla tipologia di uso del suolo e usate per stimare i flussi di deposizione provengono direttamente da CALMET.

LAPMOD può interpolare i campi meteorologici orari di CALMET con frequenza specificata in input dall'utente.

Movimento delle particelle all'interno dello strato di rimescolamento

Ogni particella si muove nel fluido atmosferico con la stessa velocità del fluido (i.e. vento medio + fluttuazione turbolenta del vento). La velocità media del vento proviene dal modello meteorologico CALMET, mentre le fluttuazioni turbolente della velocità del vento vengono calcolate da LAPMOD mediante l'equazione di Langevin con coefficienti che dipendono dalle condizioni locali di stabilità attraverso le variabili di scala dello strato limite planetario ( `u_star, w_star, L, H`).

`d u_i = a_i(x,u'_i,t) dt + b_(ij) (x,u'_i,t) d xi_j(t) `

`u_i ` rappresenta la componente media della velocità del vento nella direzione i-esima e `u'_i` rappresenta la fluttuazione turbolenta della velocità nella direzione i-esima.

La forzante casuale `d xi` viene estratta da una distribuzione normale con deviazione standard pari a `dt`.

La traiettoria di ciascuna particella viene ricostruita mediante una successione di diversi spostamenti.

`x_i (t+Delta t) = x_i (t) + Delta t (u_i + u'_i)`

Si tratta di un processo di Markov del primo ordine; ciò significa che in un flusso medio stazionario, la posizione della particella dipende solo dalla posizione che la particella occupava al tempo precedente.

Le tre componenti del moto non sono correlate tra loro.

In LAPMOD, le condizioni di stabilità da cui dipendono i coefficienti a e b dell'equazione di Langevin vengono classificate in base al numero di Richardson `L/z_i` (dove `L` è lunghezza di Monin-Obukov e `z_i` è l'altezza dello strato di rimescolamento).

Il coefficiente a per la componente verticale

Il coefficiente `a` per la componente verticale è dato da (Franzese et al., 1999):

`a = alpha w^2 + beta w + gamma`

Indicando con `W_n` il momento di ordine n della distribuzione delle fluttuazioni di velocità verticale (che dipendono da z) e con `GW_n` il corrispondente gradiente verticale, i coefficienti `alpha`, `beta` e `gamma` per condizioni instabili (`L/z_i<-1`) sono:

`alpha(z) = ( 1/3 GW_4 - 1/2 (GW_3) / W_2 (GW_3 - C_0 epsilon) - W_2 GW_2 ) / ( W_4 - W_3^2 / W_2 - W_2^2 ) `

`beta(z) = ( GW_3 - 2 alpha W_3 - C_0 epsilon) / ( 2 W_2 ) `

`gamma(z) = GW_2 - alpha W_2 `

`C_0` è la costante di Kolmogorov, pari a 3 secondo Du (1997) e `epsilon` è il rateo di dissipazione dell'energia turbolenta.

Seguendo Franzese et al. (1999), il momento secondo e il momento terzo della distribuzione sono descritti dalle seguenti espressioni analitiche:

`W_2 / w_star^2 = a_1 + a_2 (z/z_i)^(2/3) (1- z/z_i)^(4/3)`

`W_3 / w_star^3 = a_3 (z/z_i)(1- z/z_i)^2`

dove `a_1` = 0.05, `a_2` = 1.7 e `a_3` = 1.1

La kurtosi è pari a 3.5 (Hibberd and Sawford, 1994), quindi:

`W_4 = 3.5 (W_2)^2`

e il rateo di dissipazione dell'energia cinetica turbolenta è dato da (Luhar and Britter, 1989; Weil, 1990):

`epsilon = 0.4 w_star^3/z_i`

Le derivate rispetto a z dei momenti che compaiono nei coefficienti `alpha`, `beta` e `gamma` sono ottenute eseguendo le analiticamente derivate di queste parametrizzazioni.

Per condizioni neutrali e stabili (`L/z_i ge -1`) i coefficienti `alpha`, `beta` e `gamma` sono:

`alpha(z) = (GW_2)^2 / (2 (W_2)^2) = (GW_2) / W_2 `

`beta(z) = - 1 / T_L `

`gamma(z) = W_2 * GW_2 `

Dove `T_L ` è il tempo di scala Lagrangiano. Anche in questo caso le derivate rispetto a z si ottengono dalle parametrizzazioni precedenti.

Il coefficiente a per la componente orizzontale

Il coefficiente `a` per la componente orizzontale in qualsiasi condizione di stabilità è dato da:

`a = - (u'_i) / T_L `

Il coefficiente b per la componente verticale

Il coefficiente `b` per la componente verticale in condizioni instabili è dato da:

`b = sqrt (C_0 epsilon) `

mentre in condizioni neutrali e stabili è dato da:

`b = W_2 sqrt (2/ T_L) `

Il coefficiente b per la componente orizzontale

Il coefficiente `b` per la componente orizzontale in qualsiasi condizione di stabilità è dato da:

`b = W_2 sqrt (2/ T_L) `

Parametrizzazione della turbolenza e tempi di scala Lagrangiani

I tempi di scala Lagrangiani e le parametrizzazioni della turbolenza dipendono dalle condizioni di stabilità. Come già detto, l'ipotesi è che la distribuzione delle fluttuazioni della velocità delle particelle sia quella delle fluttuazioni della velocità del vento e che quindi i momenti delle distribuzioni siano gli stessi. Per il momento secondo quindi `W_2 = sigma^2` .

Tutte le equazioni presentate in questo paragrafo seguono la formulazione di Hanna et al. (1982), con l'eccezione della formulazione di Hurley and Physik (1993) per il tempo Lagrangiano nella direzione verticale in condizioni convettive.

Condizioni instabili

In condizioni instabili, i tempi di scala Lagrangiani e le fluttuazioni di velocità sono calcolati come:

`T_(Lu) = T_(Lv) = 0.15 z_i/ sigma_v`

`T_(Lw) = 0.6 z_i/ w_star`

`sigma_v = u_star(12-0.5z_i/L)^(1/3) `

Condizioni neutrali

In condizioni neutrali, i tempi di scala Lagrangiani e le fluttuazioni di velocità sono calcolati come:

`sigma_u = 2 u_star exp (- (3 f z)/(u_star)) `

`sigma_v = sigma_w = 1.3 u_star exp (- (2 f z)/(u_star)) `

`T_(Lu) = T_(Lv) = T_(Lw) = (0.5 z) / (sigma_w (1 +15 (f z)/(u_star))) `

Dove `f` è il parametro di Coriolis, z è l'altezza sopra il suolo e `u_star` è la velocità di frizione.

Condizioni stabili

In condizioni stabili, i tempi di scala Lagrangiani e le fluttuazioni della velocità sono dati da:

`sigma_u = 2 u_star ( 1 - z / z_i) `

`sigma_v = sigma_w = 1.3 u_star ( 1 - z / z_i) `

`T_(Lu) = 0.15 (z_i / sigma_u) sqrt ( z / z_i) `

`T_(Lv) = 0.07 (z_i / sigma_u) sqrt ( z / z_i) `

`T_(Lv) = 0.10 (z_i / sigma_u) sqrt ( z / z_i) `

Movimento delle particelle sopra lo strato di rimescolamento

Una particella può raggiungere e superare l’altezza dello strato di rimescolamento, sia a causa del suo moto, sia perché lo strato di rimescolamento si abbassa al variare delle condizioni meteorologiche.

In questi casi, se una particella era all’interno dello strato di rimescolamento durante l’intervallo di tempo precedente, la sua fluttuazione di velocità verticale è posta pari a zero; altrimenti, non vengono effettuate modifiche.

Sopra lo strato di rimescolamento, le fluttuazioni di velocità delle particelle seguono per le tre componenti l’equazione di Langevin per condizioni omogenee:

`a = - 1 / T_L `

`b = W_2 sqrt (2/ T_L) `

dove `T_L` è posto pari a 1000 s e `W_2` pari a 0.01 m2s-2. Si assume assenza di correlazione tra le componenti.

Riflessione della particelle al suolo

Le particelle, quando arrivano al suolo (cioè quando sono al di sotto della lunghezza di rugosità `Z_0`), subiscono una riflessione verticale perfetta, in base alla relazione:

` Z_(n ew) = 2 Z_0 - Z_(old) `

Anche lo spostamento dovuto alla fluttuazione di velocità verticale subisce una riflessione perfetta in maniera del tutto analoga.

Effetti dovuti al galleggiamento termico

I camini possono essere considerati come sorgenti puntuali con un effetto di galleggiamento termico. Questo effetto può essere stimato seguendo l’approccio proposto da Janicke e Janicke (2001) oppure da Webster e Thomson (2005).

Algoritmo di Janicke e Janicke

Uno degli algoritmi contenuti in LAPMOD per la descrizione del plume rise è quello descritto da Janicke e Janicke (2001), che ha la peculiarità di essere adatto anche per piume con elevato contenuto di umidità.

L’algoritmo è tridimensionale e può quindi essere applicato anche a situazioni con campo anemologico molto complesso, e non richiede semplificazioni poco realistiche come velocità del vento dipendente solo dalla quota verticale e direzione costante, o addirittura vento costante per tutta l'altezza della piuma. Inoltre il modello non adotta l'ipotesi di Boussinesq (cioè non assume che la differenza tra la densità della piuma e la densità dell'atmosfera sia molto minore della densità dell'atmosfera), quindi può essere applicato anche a situazioni in cui la temperatura di emissione è molto maggiore della temperatura atmosferica.

Ciò significa che l’algoritmo è adatto per una vasta gamma di tipologie di sorgenti, dalle torri di raffreddamento, dove la temperatura di emissione differisce da quella atmosferica solo per qualche decina di gradi, ai camini dotati di scrubber (e.g. desolforatori) dove la temperatura della piuma può differire da quella atmosferica anche per centinaia di gradi.

La formulazione matematica del modello permette inoltre di simulare un camino con una direzione arbitraria di uscita e di considerare sia piume umide (sature o non sature) sia piume secche.

Di seguito vengono riportate le equazioni principali del modello di Janicke e Janicke (2001). L’algoritmo viene risolto numericamente all’interno di LAPMOD integrando le equazioni per mezzo di un metodo Runge-Kutta al quarto ordine.

Tutte le equazioni del modello vengono scritte in funzione della distanza lungo l'asse della piuma, denotata con s. Le relazioni tra le coordinate cartesiane e la coordinata s sono le seguenti:

` dx / (ds) = u_x / u; dy /(ds) = u_y / u; dz/(ds) = u_z/u`

dove `u` la velocità della piuma misurata lungo il suo asse e `u_x`,`u_y` e `u_z` sono le componenti cartesiane della velocità della piuma.

La conservazione della massa d'aria secca è descritta dall'equazione

` (dM) / (ds) = epsilon`

Il flusso di massa M è dato da:

` M = A rho u`

dove A è la sezione della piuma e `rho` è la densità della piuma. Il coefficiente di entrainment `epsilon`, come verrà indicato nel seguito, può essere descritto attraverso tre diverse formulazioni all’interno di LAPMOD.

La conservazione del momento è descritta attraverso l'equazione

` (d(M u_i)) / (ds) = epsilon + delta_(ij) A g (rho_a -rho)`

dove il pedice "i" indica una delle 3 direzioni cartesiane (i=1 x, i=2 y, i=3 z) e ` delta_(ij) ` è il delta di Kronecker, quindi il secondo addendo a destra dell'uguale è diverso da zero solo per la componente verticale. Il pedice “a” indica una proprietà dell'atmosfera.

Il contenuto d'acqua in fase vapore della piuma viene indicato attraverso l'umidità specifica q, definita come il rapporto tra la massa di vapore acqueo e la massa totale di aria umida in uno specifico volume. Il contenuto di acqua liquida viene indicato con `eta` ed è il rapporto tra la massa d'acqua in fase liquida e la massa totale di aria umida in uno specifico volume. Il contenuto totale di acqua nella piuma (fase vapore più fase liquida) è quindi dato da `xi = eta + q `. L'equazione per la conservazione della massa d'acqua è:

` (d(M xi)) / (ds) = epsilon + xi_a`

Il passaggio di stato della massa d'acqua in un punto della piuma viene simulato calcolando il valore di saturazione dell'umidità specifica in funzione della temperatura in quel punto della piuma. Se il contenuto di vapore specifico è superiore al valore di saturazione, la parte in eccedenza condensa aumentando la massa in fase liquida.

L'equazione di bilancio per l'entalpia è data da:

` d/(ds) (M(T- h_v/c_p eta)) = -M g/c_p u_z /u rho_a / rho + epsilon (T_a - h_v/c_p eta_a)`

dove `h_v` è il calore latente di evaporazione dell'acqua è `c_p` è la capacità termica specifica a pressione costante.

Le equazioni descritte vengono risolte numericamente dopo aver specificato opportune condizioni iniziali. Le soluzioni forniscono il valore medio delle variabili di interesse su un segmento ortogonale all'asse della piuma e per un certo valore della variabile s. Si assume che le variabili scalari (temperatura, contenuto d'acqua e concentrazione di eventuali inquinanti) abbiano un profilo ortogonale all'asse della piuma della forma:

`X = X_a + X^star exp (- r^2 / b^2) `

dove X è la generica variabile scalare, r è la distanza dall'asse della piuma e le quantità `X^star` e b sono date da soddisfare:

`X^star = 3 (X - X_a) `

`b = R / sqrt 3 `

Una applicazione dell’algoritmo al problema della visibilità di piume umide è stata descritta in Presotto et al. (2005). Infatti, una volta calcolata la temperatura media della piuma alla distanza assiale s, il profilo di temperatura si ottiene con una equazione come quella indicata sopra, ed è perciò possibile calcolare l'umidità specifica di saturazione e la regione dove l'umidità relativa della piuma è pari al 100%, cioè dove avviene la condensazione. La regione della piuma dove il contenuto specifico di acqua liquida `eta` supera un dato valore di soglia costituisce la porzione visibile della piuma.

Per quanto riguarda il coefficiente di entrainment ε, LAPMOD permette di scegliere tra tre diverse formulazioni: Janicke and Janicke (2001), Rezacova e Sokol (2000) e Orville et al. (1980).

Il coefficiente di entrainment proposto da Janicke and Janicke (2001) è:

`epsilon = 2 pi rho_a E `

con `rho_a` densità dell'atmosfera (anche nel seguito il pedice "a" indicherà una proprietà dell'atmosfera), ed E funzione di entrainment descritta come:

`E = lambda R (alpha/2 (w_(per))^2 / w + beta (w_(par))^2 / w + u gamma /(Fr^2))`

dove R è il raggio della piuma, `lambda`, `alpha`,`beta`,`gamma`, sono coefficienti empirici, w è il modulo della velocità relativa della piuma (velocità del vento meno velocità della piuma), `w_(par)` e `w_(per)` sono le componenti della velocità relativa parallela e perpendicolare all'asse della piuma, `Fr^2` è il quadrato del numero di Froude calcolato come:

`Fr^2 = (rho u^2) / (abs (rho_a - rho) g R)`

Il coefficiente di entrainment proposto da Rezacova e Sokol (2000) è dato da

`epsilon = (alpha_e T_v) / (R T_(Va)) M`

dove `alpha_e` è un coefficiente empirico (`alpha_e` = 0.3), R è il raggio della piuma (m), `T_V` è la temperatura virtuale della piuma (K), `T_(Va)` è la temperatura virtuale dell’atmosfera (K) e M è il flusso di massa (kg m‑1 s‑1). Sostituendo l’espressione del flusso di massa nell’equazione precedente si ha:

`epsilon = pi alpha_e rho u T_v / T_(Va) R`

Orville et al. (1980) adattano un modello per la formazione delle nubi alla descrizione della piuma formata dalle emissioni di torri di raffreddamento. Il coefficiente di entrainment per una massa d’aria che sale verticale viene descritto da:

`epsilon = (2 alpha) / R M`

con `alpha = 0.075` coefficiente empirico. M nell’equazione precedente si ottiene

`epsilon = 2 alpha pi rho u R`

Algoritmo di Webster e Thomson

Di seguito vengono riportate le equazioni principali del modello di Webster e Thomson (2002). L’algoritmo viene risolto numericamente all’interno di LAPMOD integrando le equazioni per mezzo di un metodo Runge-Kutta al quarto ordine.

Tutte le equazioni del modello vengono scritte in funzione della distanza lungo l'asse della piuma, denotata con s. Le relazioni tra le coordinate cartesiane e la coordinata s sono le seguenti:

` dx / (ds) = u_x / u; dy /(ds) = u_y / u; dz/(ds) = u_z/u`

dove `u` la velocità della piuma misurata lungo il suo asse e `u_x`,`u_y` e `u_z` sono le componenti cartesiane della velocità della piuma.

Le equazioni di conservazione accoppiate risolte dal modello sono le seguenti:

` (dF_(m))/(ds) = epsilon`

` (dF_(Mx))/(ds) = -F_m (du_a)/(ds) - D_x`

` (dF_(My))/(ds) = -F_m (dv_a)/(ds) - D_y`

` (dF_(Mz))/(ds) = -F_m (dw_a)/(ds) +B - D_z`

` (dF_H)/(ds) = -F_m c_p (d theta_a)/(ds)`

dove ` ul D = (D_x,D_y,D_z)` è la forza di drag, B è la forza di buoyancy, ` epsilon` è il rateo di entrainment della massa d'aria ambiente, ` ul (u_a) = (u_a,v_a,w_a)` è la velocità media del vento, ` c_p` è il calore specifico dell'aria a pressione costante. Il flusso di massa ` F_m` , il flusso di momento ` ul (F_M) = (F_(Mx),F_(My),F_(Mz))` e il flusso di calore ` F_H` sono dati da:

` F_m = pi b^2 rho_p u_s`

` F_(Mi) = (u_(p i) - u_(ai)) F_m`

` F_H = c_p (theta_p - theta_a) F_m`

dove ` b=b(t)` è il raggio della piuma, ` u_s` è il modulo della velocità della piuma ` ul u_p = (u_p, v_p, w_p)` . I pedici p ed a denotano una proprietà della piuma e dell'aria rispettivamente, mentre il pedice i denota una delle tre componenti cartesiane.

La forza di galleggiamento B e il rateo di entrainment ` epsilon` vengono calcolati come:

` B = pi b^2 g (rho_a - rho_p)`

` epsilon = 2 pi b u_e`

dove ` u_e` è la velocità di entrainment ottenuta come:

` u_e = u_e^(rise) + u_e^(turb)`

` u_e^(rise) = alpha_1 abs (Delta u_s) + alpha_2 abs (Delta u_N)`

` u_e^(turb) = alpha_3 min ((epsilon b)^(1/3),alpha_w (1+t/(2 T_L))^(-1/2))`

dove, la prima componente ` u_e^(rise)` è dovuta al moto relativo tra la piuma e l'ambiente, la seconda componente `u_e^(turb)` è dovuta alla turbolenza ambientale, `alpha_1`, `alpha_2` e `alpha_3` sono costanti, `alpha_w` delle fluttuazioni verticali di velocità, `Delta u_s` e ` Delta u_N` sono le componenti delle velocità relative nelle direzioni parallela e ortogonale all'asse della piuma, ` T_L` è la scala di tempo Lagrangiano.

infine, le tre componenti cartesiane della forza di drag sono date ciascuna da:

` D_i = pi b rho Delta u_s abs (Delta u_N) c_D `

dove `c_D ` è il coefficiente di drag.

I valori utilizzati per le costanti dell'algoritmo di plume rise sono tipicamente `alpha_1 = 0.11`, `alpha_2 = 0.5`, `alpha_3 = 0.655` e `c_D = 0.21`. Tali valori possono tuttavia variare, e sono quindi input del modello LAPMOD. I coefficienti di entrainment dovuti al moto relativo tra piuma e ambiente utilizzati in ADMS sono ad esempio `alpha_1 = 0.057` e `alpha_2 = 0.50`. Il valore di `alpha_1` è molto diverso dai valori utilizzati in altri modelli, e in condizioni di calma di vento non permette di determinare un'altezza di plume rise paragonabile a quella delle formule classiche di Briggs o Weil. Il modello TAPM invece utilizza `alpha_1 = 0.100 ` e `alpha_2 = 0.600`. Valori tipici dei due coefficienti di entrainment per moto relativo della piuma sono riportati in Tabella 2 (da Webster e Thomson, 2002). Per il coefficiente di entrainment dovuto alla turbolenza atmosferica `alpha_3` e per il coefficiente di drag `c_D` non esistono attualmente molti valori in letteratura, quindi si utilizzano i valori indicati in precedenza.

Condizioni meteorologiche	Prevalenza del momento	Prevelenza del galleggiamento
Forte vento ambientale (Bent over plume). Nessuna stratificazione.	`alpha_2 = 0.35` (Briggs) `alpha_2 = 0.6` (Weil)	`alpha_2 = 0.61` (Briggs) `alpha_2 = 0.6` (Weil)
Nessun vento ambientale (piuma verticale). Nessuna stratificazione.	`alpha_1 = 0.11` (Briggs)	`alpha_1 = 0.125` (Briggs) `alpha_1 = 0.11` (Weil)
Forte vento ambientale (Bent over plume). Stratificazione stabile.	`alpha_2 = 0.6` (Briggs)	`alpha_2 = 0.6` (Briggs)

Terminali dei camini

La geometria del terminale determina direzione e quantità di moto iniziali del pennacchio. Come descritto da Bellasio e Bianconi (2025), uno sbocco verticale libero fornisce una spinta iniziale verso l’alto; il galleggiamento termico dipende invece dalla differenza di temperatura tra fumi e ambiente. Gli sbocchi orizzontali o coperti sopprimono la quantità di moto verticale iniziale, rendendo particolarmente importante il galleggiamento termico.

Per sbocchi inclinati e orizzontali, azimut e inclinazione definiscono la direzione del rilascio. La velocità relativa tra pennacchio e vento controlla l’entrainment: un maggiore mescolamento porta più rapidamente le proprietà del pennacchio verso quelle ambientali. Gli sbocchi orizzontali mantengono quindi il loro orientamento effettivo, senza assumere un’emissione sottovento.

Per una copertura antipioggia, la formulazione dell’articolo raddoppia il diametro iniziale, assegna una velocità verticale di 0.001 m/s e devia il flusso orizzontalmente lungo il vento a un quarto della velocità di uscita dal camino. Sotto 0.5 m/s di vento, la direzione orizzontale viene scelta casualmente. Queste ipotesi rappresentano l’espansione e la deviazione sotto la copertura.

Il terminale influenza altezza del pennacchio, diluizione e concentrazioni vicine. Lo studio ha trovato gli impatti minori ai recettori più prossimi per sbocchi verticali liberi, soprattutto con basse temperature dei fumi; le differenze dipendono anche da meteorologia e orografia.

Turbolenza indotta dal plume rise

L'effetto del plume rise è quello di innalzare l'inquinante e di aumentare il livello di turbolenza all'interno della piuma. Il livello di turbolenza all'interno della piuma è elevato in prossimità del rilascio, dove temperatura e velocità di emissione sono elevate, e diminuisce allontanandosi dal punto di emissione a causa dell'entrainment di aria con caratteristiche dell'ambiente circostante.

Come indicato in precedenza, una particella computazionale in atmosfera si muove per effetto del vento medio e di fluttuazioni turbolente. Una particella all'interno di una piuma in fase di risalita viene invece spostata utilizzando, al posto delle componenti del vento medio, le componenti della velocità della piuma, calcolate da LAPMOD per mezzo degli algoritmi di Janicke e Janicke o di Webster e Thomson, a seconda delle indicazioni dell'utente. A tali componenti vengono in ogni caso sommate le fluttuazioni turbolente. Per considerare la turbolenza indotta dal plume rise, alla posizione della particella viene inoltre sommato uno spostamento casuale in ogni direzione. L'ampiezza dello spostamento casuale viene determinata estraendo un valore da una distribuzione Gaussiana di media nulla e varianza V data da (Webster e Thomson, 2002):

`V = (r^2 (t+Delta_t) - r^2(t)) * 0.25`

dove r è il raggio della piuma, t indica il tempo attuale e `t+Delta_t` il tempo futuro. L'ampiezza dello spostamento è quindi elevata in prossimità del rilascio, quando l'attività di entrainment è massima e il raggio della piuma cresce molto rapidamente, mentre si riduce allontanandosi dalla sorgente, sino ad annullarsi al termine del plume rise. L'ampiezza dello spostamento diminuisce quindi con il livello di turbolenza indotto dalla piuma. Al termine del plume rise la particella torna a muoversi per effetto del vento medio e delle fluttuazioni turbolente dovute ai soli effetti atmosferici.

Stack tip downwash

Per le sorgenti puntiformi con galleggiamento termico e/o meccanico, LAPMOD considera anche lo stack tip downwash (STD). Il fenomeno dello STD può avvenire quando il rapporto tra la velocità di uscita dei fumi e la velocità del vento alla bocca del camino è piccolo (cioè quando la velocità del vento è maggiore della velocità di uscita dei fumi). In tal caso la piuma può essere catturata nella scia del camino, con il risultato di aumentare i valori di concentrazione immediatamente sottovento al camino. Tale fenomeno viene simulato per mezzo di un algoritmo di Briggs che ha l’effetto pratico di diminuire l’altezza fisica del camino quando le condizioni sono tali da generare STD. Quando il rapporto tra la velocità di uscita dei fumi (`v_(out)`) e la velocità del vento alla bocca del camino (`w`) è minore di 1.5 l’altezza del camino viene calcolata come:

`H'_S = H_S + 2 D ( v_(out) / w - 1.5)`

dove `H'_S` è l’altezza reale del camino e `D` è il suo diametro. Se invece il rapporto tra le velocità è maggiore o uguale a 1.5, l’altezza del camino risulta invariata (cioè non si osserva lo STD).

Partial plume penetration

Le figure nella tabella seguente mostrano schematicamente come può comportarsi la piuma emessa da un camino relativamente all’altezza di rimescolamento Hmix. Si possono presentare quattro situazioni.

	Piuma emessa direttamente sopra Hmix. Avviene tipicamente per camini molto elevati e in condizioni di stabilità atmosferica, quando l’altezza di rimescolamento è abbastanza bassa, tipicamente in condizioni notturne. In queste condizioni la bocca del camino si trova sopra Hmix, e la piuma viene emessa all’interno di uno strato caratterizzato da scarsa turbolenza, quindi è in grado di muoversi compatta fino a distanze relativamente grandi.
	Piuma emessa sotto Hmix con sufficiente galleggiamento da “bucare” Hmix. Avviene per altezze di rimescolamento non molto elevate e per camini con valori elevati del parametro di galleggiamento, caratterizzati cioè da grandi diametri di emissione, elevate velocità di uscita e temperature di emissione molto elevate rispetto alla temperatura ambiente.
	Piuma intrappolata sotto Hmix. Avviene nelle condizioni in cui il parametro di galleggiamento non è molto elevato, mentre sono elevate la velocità del vento alla quota di emissione del camino, la differenza di altezza tra lo strato di rimescolamento e il camino, e la forza dell’inversione termica in quota (cioè la differenza di temperatura attraverso Hmix).
	Piuma parzialmente intrappolata sotto Hmix. Avviene nelle condizioni in cui il parametro di galleggiamento è confrontabile con l’effetto combinato della velocità del vento alla quota di emissione del camino, della differenza di altezza tra lo strato di rimescolamento e il camino, e della forza dell’inversione termica in quota.

LAPMOD è in grado di simulare le quattro condizioni descritte. In particolare l’intrappolamento della piuma dovuto alla presenza di una inversione in quota viene simulato per mezzo dell’algoritmo descritto da Manins (1979), lo stesso utilizzato anche nel modello CALPUFF. Viene definito un parametro di penetrazione P come:

` P = (F_b) / (u_s b_i (H_(mix) - H_s))`

Dove ` u_s` è la velocità del vento (m/s) alla quota della bocca del camino, ` H_s` è l’altezza del camino (m), ` H_(mix)` è l’altezza di rimescolamento (m), ` b_i` è la forza dell’inversione (m/s2) e ` F_b` è il parametro di galleggiamento (m4/s3). La forza dell’inversione è data da:

` b_i = DeltaT / T_(as)`

Dove ` DeltaT` è la differenza di temperatura attraverso l’inversione (K), ` T_(as)` è la temperatura atmosferica alla quota della bocca del camino (K), e ` g` è l’accelerazione di gravità (m/s2). Il parametro di galleggiamento ` F_b` è dato da:

` F_b = (g d^2 v_(out) (T_(out) - T_(as)))/ (4 T_(as))`

dove ` d` è il diametro del camino (m), ` v_(out)` è la velocità di emissione (m/s) e ` T_(out)` è la temperatura di emissione (K).

Noto il parametro di penetrazione ` P` , si calcola la frazione ` f` della massa della piuma che rimane sotto ` H_(mix)` . Se ` P lt 0.08` , ` f = 1` , cioè la piuma rimane completamente sotto lo strato di inversione, mentre se ` P gt 0.3` , ` f = 0` , cioè la piuma buca completamente lo strato di inversione in quota. Negli altri casi, cioè quando ` 0.08 lt P le 0.3` , la frazione della massa della piuma che rimane sotto Hmix viene calcolata come:

` f = 0.08 / P - P + 0.08`

L’andamento della frazione di massa intrappolata al variare del parametro di penetrazione viene riportato nel grafico seguente.

Nel caso in cui ` f > 0` , la quota finale di plume rise viene limitata al valore:

` z_p = H_s + (1- f/3) (H_(mix) - H_s)`

Il parametro di penetrazione (` P` ), la frazione di massa che rimane sotto Hmix (` f` ) e l’eventuale valore limite del plume rise (` z_p`) vengono calcolati all’interno della routine PAPLPE. Nel caso in cui da tale routine risulti una frazione ` f > 0` , le limitazioni alla buoyancy e al plume rise delle particelle vengono imposte all’interno della routine PMOVE.

Nel caso in cui ` 0 lt f lt 1` , cioè di nube solo parzialmente intrappolata, per ogni particella emessa dal camino viene generato un numero casuale (` R` ) compreso tra 0 e 1, e la particella viene bloccata se ` R le f` , mentre è libera di salire oltre ` H_(mix)` nel caso contrario.

Riepilogando, nella routine PAPLPE vengono gestiti i seguenti casi:

` Hs gt Hmix`	Tutta la piuma è emessa fuori dallo strato limite.
` DeltaT = T(z_2) - T(z_1) lt 0`	Non è inversione, la routine non è applicabile, la piuma segue esattamente il plume rise senza modifiche.
` Fb lt 0 `	Rilascio freddo (a temperatura minore della temperatura ambiente), la routine non è applicabile, la piuma segue esattamente il plume rise senza modifiche.
` P gt 0`	La piuma buca completamente Hmix.
` P lt 0.08 `	La piuma rimane completamente sotto Hmix.
` 0.08 le P lt 0.3`	La piuma buca parzialmente Hmix. La frazione di massa che rimane sotto Hmix è data da f.

Building downwash

L?implementazione del building downwash in LAPMOD ? descritta da Bellasio e Bianconi (2026).

Edifici e ostacoli modificano il flusso atmosferico circostante. Un pennacchio che entra nella regione perturbata può raggiungere il suolo più rapidamente, aumentando le concentrazioni vicino alla sorgente. Il building downwash è quindi rilevante per gli odori e i recettori vicini. È distinto dallo stack-tip downwash, causato dal flusso intorno al camino stesso.

Si distinguono due regioni (Schulman et al., 2000). La cavità, o scia vicina, è delimitata dalle linee di corrente che si separano dal tetto e dai lati dell’edificio e dal loro successivo riattacco. Presenta basse velocità medie e moto intenso, spesso ricircolante. Oltre la cavità si estende la scia lontana: il vento medio resta ridotto e la turbolenza aumentata, con disturbi progressivamente più deboli verso le condizioni ambientali.

Le dimensioni dipendono dall’altezza H, dalla larghezza W proiettata trasversalmente al vento e dalla lunghezza L lungo il vento. La scala caratteristica di Wilson (1979) descrive la geometria, compresi riattacco sul tetto ed estensione sottovento. Poiché le dimensioni proiettate cambiano con la direzione del vento, anche la regione interessata dipende dall’orientamento del flusso rispetto all’edificio.

LAPMOD usa le relazioni empiriche di Petersen e Guerra (2018), ricavate da prove in galleria del vento per strutture rettangolari e aerodinamiche, come serbatoi e torri di raffreddamento a tiraggio naturale. Le relazioni distinguono gli incrementi laterale e verticale della turbolenza, generalmente minore quello laterale, e descrivono il rapido ritorno alle condizioni ambientali sopra il tetto. Le strutture aerodinamiche producono scie e deficit di velocità inferiori ai blocchi rettangolari. Una maggiore rugosità a monte riduce l’intensità di turbolenza aggiuntiva associata alla scia. Questi sviluppi sono parte delle ricerche su PRIME2 (Petersen e Guerra, 2018; Petersen et al., 2022).

La scia modifica trasporto medio e moto turbolento delle particelle nella regione interessata. Il trattamento si applica a sorgenti puntiformi, areali e volumetriche, anche dopo la fine del plume rise. L’intensità varia nello spazio: gli effetti al suolo dipendono da posizione e quota della sorgente e geometria degli edifici. L’uso di queste relazioni non implica l’adozione dell’intero schema di plume rise PRIME2: LAPMOD mantiene le formulazioni di Janicke e Janicke (2001) e Webster e Thomson (2002).

Vedere BUILDINGS.TXT per l’input e gli aggiornamenti dell’implementazione per il percorso PRIME2 e la correzione della direzione del vento.

Intervalli di integrazione temporale

Ciascuna particella viene spostata indipendentemente e la sua posizione viene aggiornata fino al momento in cui si procede al calcolo della concentrazione e della deposizione. Seguendo Thomson (1987), la condizione da applicare per determinare gli intervalli di integrazione brevi per spostare ogni particella lungo la verticale è:

` Deltat = min(0.05 T_L, 0.1 W_2 / abs(a). 0.01 W_2 / abs (w' GW_2))`

Per le componenti orizzontali la condizione sull’intervallo di integrazione è ` DeltaT = 0.15 T_L` , come proposto dall’analisi di Wilson e Zhuang (1989) sugli errori dovuti alla discretizzazione con differenti intervalli di integrazione in condizioni omogenee.

Inoltre, viene applicata la condizione di Courant, così che ogni particella risente delle condizioni meteorologiche di tutte le celle che ha attraversato nel suo percorso.

Decadimento radioattivo

Ad ogni intervallo di integrazione ` Deltat` , la massa di ciascuna particella in atmosfera è diminuita in base al coefficiente di decadimento radioattivo ` lambda ` del radionuclide:

` M(t+Deltat) = M(t) exp (-lambda Deltat)`

Il decadimento radioattivo viene applicato anche alla massa di radionuclidi depositati al suolo.

Non vengono considerate le catene di decadimento e non vengono generate nuove specie e/o nuove particelle.

Deposizione

Una particella all’interno dello strato superficiale (stimato come il 10% dell’estensione dello strato di rimescolamento) può perdere una parte della sua massa a causa della deposizione secca o umida. Questo avviene sia per i gas che per gli aerosol.

In caso di deposizione secca, la perdita di massa in un dato intervallo di tempo ` Deltat ` dipende esponenzialmente dalla velocità di deposizione, ` v_d`:

` M(t+Deltat) = M(t) exp (-(v_d Deltat)/Z_S)`

dove ` Z_S` è l’altezza dello strato superficiale.

La velocità di deposizione per i gas ` v_d` si ricava definendo delle resistenze, in analogia con i circuiti elettrici:

` v_d = 1/ (r_a + r_b + r_c) `

La resistenza aerodinamica ` r_a` viene calcolata come suggerito da Seinfeld e Pandis (1998):

Condizioni di stabilità
Atmosfera stabile	` r_a = 1/(k u_(star)) (log(z/z_0) + 4.7 (xi - xi_0))`
Atmosfera neutrale	` r_a = 1/(k u_(star) (log(z/z_0))`
Atmosfera instabile	` r_a = 1/(k u_(star)) (log(z/z_0) + log ( ((eta_0^2 +1)(eta_0 + 1)^2) / ((eta_r^2 +1)(eta_r + 1)^2) )) + 2(tan^(-1) eta_r - tan^(-1)i eta_0 ) )`

dove ` eta_0 = (1-15 z_0/L)^(1/4)`,` eta_r = (1-15 z_r/L)^(1/4)`, `z_0` è la lunghezza di rugosità e `z_r` è l'altezza di riferimento.

La resistenza superficiale è calcolata seguendo Maryon e Ryall (1996):

` r_b = 8 / u_(star)`

La resistenza del canopy superficiale, ` r_c ` è assunta pari a 500 s/m (Verver e De Leeuw, 1992, Baklanov e Soerensen, 2001).

La velocità di deposizione secca ` v_d` per gli aerosols è determinata considerando la velocità di sedimentazione ` v_(Se t tl)` e gli algoritmi descritti da Zhang e al. (2001). La velocità di sedimentazione è calcolata in base alla legge di Stokes:

` v_(Se t t l) = (rho d^2 g) / (18 eta) C(d) `

dove ` r` è la densità di aerosol (g/m3), ` g` l’accelerazione di gravità (m/s2), ` eta ` è la viscosità dell’aria (g/(m s)) e ` C(d)` è il fattore di correzione di Cunningham che dipende dal diametro dell’aerosol (determinato come spiegato in Seinfeld e Pandis, 1998).

La velocità di deposizione dell’aerosol ` v_d` è calcolata dall’algoritmo di Zhang e al. (2001), che aumenta la velocità di deposizione ` v_(Se t t l)` di un fattore calcolato considerando la categoria di occupazione del suolo, la velocità di frizione ` u_(star)` e la stagione dell’anno. Le categorie di occupazione del suolo considerate da questo algoritmo sono sintetizzate neòlla tabella seguente:

ID	Descrizione
1	Urban
2	Agricultural
3	Range land
4	Deciduous forest
5	Coniferous forest
6	Mixed forest including wetland
7	Water, both salt and fresh
8	Barren land, mostly desert
9	Non forested wetland
10	Mixed agricultural and range land
11	Rocky open areas with low growing shrubs

La velocità di deposizione secca per gas e aerosol è calcolata in LAPMOD dalla routine DEPOSD e dalle sue subroutines.

L’attuale versione di LAPMOD determina la deposizione umida solo per gli inquinanti in fase di aerosol. Questa può essere considerata una buona approssimazione, in quanto la deposizione umida è molto più efficace per gli aerosol che per i gas. La velocità di deposizione umida è determinata come la somma del dilavamento dovuto al rainout e a quello dovuto al washout. L’efficienza delle gocce d’acqua nel catturare gli aerosol dipende da molti fattori:

• l’impatto dell’aerosol sulla goccia d’acqua,
• l’intercettazione delle particelle da parte della goccia d’acqua,
• il moto Browniano delle particelle nella goccia d’acqua,
• il processo di enucleazione di una goccia d’acqua da parte di una particella,
• l’attrazione elettrica,
• l’attrazione termica
• la diffusoforesi

La velocità di deposizione umida ` v_w` di un aerosol si determina come descritto da Baklanov e Sorensen (2001):

` v_w = Lambda H_C`

dove ` Lambda` è il coefficiente di washout medio sulla verticale e ` H_C` è l’altezza della base della nuvola. Il coefficiente di washout è una funzione del raggio dell’aerosol e del tasso di precipitazione come descritto dalle seguenti equazioni:

` r lt 1.4 mu m`	` Lambda = a_0 q^(0.79 `
`1.4 mu m lt r lt 10 mu m`	`Lambda = (b_0+b_1 r + b_2 r^2 + b_3 r^3)(a_1 q + a2 q^2)`
`r ge 10 mu m`	`Lambda = a_1 q + a_2 q^2 `

dove `a_0`, `a_1`, `a_2`, `b_0`, `b_1`, `b_2`, e `b_3` sono costanti, `r` è il raggio dell’aerosol (`mu m`) e `q` è il tasso di precipitazione (mm/h). Per il processo di rainout il coefficiente di dilavamento è calcolato in accordo con la precedente equazione usando la formula per il raggio più piccolo.

Le velocità di deposizione umida per la fase aerosol sono calcolate in LAPMOD dalla routine DEPOSW e dalle sue subroutines.

La variazione di massa dovuta alla deposizione totale (secca e umida) in un intervallo può essere espressa dalla relazione:

`M(t+Delta_t) = M(t) exp (- nu Delta_t)`

dove ` nu = d + w` è la somma delle velocità di deposizione secca e umida, ciascuna divisa per la propria caratteristica lunghezza di scala.

Anche se la deposizione secca e umida sono processi in competizione, devono essere calcolati separatamente per valutare il loro contributo alla deposizione totale. Per questa ragione è necessario calcolare le due masse equivalenti ` M_d` e ` M_w` ciascuna sottoposta separatamente ad uno dei due processi:

`M_d(t+Delta_t) = M_(star)(t) exp (- d Delta_t) `

`M_w(t+Delta_t) = M_(star)(t) exp (- w Delta_t) `

Poichè `M_d(t+Delta_t) = M_d(t+Delta_t) + M_w(t+Delta_t) ` si ricava che la massa equivalente (`M_(star)`) che deve essere usata nelle due equazioni sopra per calcolare separatamente la deposizione secca e umida in un intervallo di tempo `Delta_t` è:

`M_(star)(t) = M(t) (2 - exp(-d Delta_t) -exp(-w Delta_t)) / (1 - exp(-(d+w) Delta_t))`

Calcolo della concentrazione

Aggiornato all’8 settembre 2026. Vedere gli aggiornamenti dei kernel, del building downwash e di OpenMP.

Molti modelli Lagrangiani a particelle determinano la concentrazione istantanea in un punto (x,y,z) con la tecnica del conteggio, cioè contando le particelle presenti in un ipotetico volume di atmosfera e dividendo la massa totale di tali particelle per il volume che le contiene:

`C_(IST)(x,y,z) = 1 / (Deltax Deltay Deltaz) sum_(i=1)^N m_i`

dove:

• `m_i ` è la massa associata alla i-esima particella computazionale relativa alla specie in esame,
• `Deltax`, `Deltay` e `Deltaz` (m) sono le tre dimensioni del volume ipotetico all'interno del quale si calcola la concentrazione,
• N è il numero di particelle computazionali presenti nel volume di atmosfera considerato.

Il metodo del conteggio è concettualmente semplice ma ha alcuni svantaggi:

• richiede la generazione di un elevato numero di particelle, quindi aumenta i tempi di simulazione;
• costringe a dividere il dominio di simulazione in celle di calcolo (come per i modelli Euleriani) e la concentrazione calcolata dipende dal volume di ogni cella, producendo campi di concentrazione con gradienti sottostimati se i volumi di conteggio sono molto grandi e forti discontinuità se i volumi sono troppo piccoli.

Un approccio alternativo consiste nell'utilizzare un Kernel density estimator.

Con il termine kernel, nel contesto dei modelli di dispersione atmosferica Lagrangiani a particelle, si indica una funzione matematica che permette di calcolare la concentrazione in ogni punto dello spazio a partire da una specifica disposizione delle particelle.

La concentrazione in corrispondenza dei recettori è calcolata come la somma dei contributi di tutte le particelle presenti all’interno del dominio delle simulazione. Tutte le particelle, durante il moto, vengono diffuse dalla turbolenza atmosferica (quindi la massa di due particelle emesse allo stesso istante, dopo un certo tempo, avrà in genere una diversa distribuzione spaziale data la natura stocastica del processo diffusivo). Per il calcolo delle concentrazioni quindi le particelle sono trattate come delle nuvole (distribuzioni) con diversa dimensione e diversa area di influenza.

Il contributo di una singola particella alla concentrazione nel punto (x,y,z) distante `d_x`, `d_y` e `d_z` dalla particella è dato da:

`C_i(x,y,z) = (eta m_i)/((2 pi)^(3/2) sigma_(x i) sigma_(yi) sigma_(zi)) K_x(d_(x i)/beta_(x i)) K_y(d_(yi)/beta_(yi)) K_z(d_(zi)/beta_(zi))`

dove `m_i` è la massa e `sigma_(x i)`, `sigma_(yi)` e `sigma_(zi)`, sono le sigma rispettivamente orizzontali e verticale della i‑esima particella e `beta_(x i)`, `beta_(yi)` e `beta_(zi)` sono delle grandezze di scala. Le funzioni `K` descrivono la distribuzione spaziale della massa associata alla particella e dal punto di vista formale devono unicamente sottostare alla condizione di integrabilità nelle tre dimensioni.

LAPMOD incorpora diverse formulazioni delle funzioni `K`.

Kernel Gaussiano (CCA1)

Il kernel gaussiano (Yamada e Bunker, 1988) utilizza una forma esponenziale per le funzioni `K`:

`C_i(x,y,z) = (eta m_i)/((2 pi)^(3/2) sigma_(hi)^2 sigma_(vi)) exp( -d_(x i)^2/(2 sigma_(hi)^2)) exp( -d_(y i)^2/(2 sigma_(hi)^2)) [exp( -d_(z i)^2/(2 sigma_(vi)^2)) + exp( -d_(z + z i)^2/(2 sigma_(zi)^2))]`

dove mi è la massa e `m_i` e `sigma_(hi)` (ipotizzando `sigma_(hi) = sigma_(x i) = sigma_(yi)` ) e `sigma_(vi)` sono le sigma rispettivamente orizzontale e verticale della i‑esima particella; `z_i` è la coordinata verticale della i-esima particella e ` z` è la coordinata verticale del recettore in esame. Nel termine verticale dell’equazione viene considerata anche la riflessione dal suolo.

Il fattore `eta` che compare nell’espressione della concentrazione garantisce la normalizzazione. Infatti, il contributo di una data particella al calcolo della concentrazione in corrispondenza di un recettore viene troncata quando la distanza supera k `sigma` dove k è un numero intero (di solito si pone k=3) e per questo motivo la distribuzione di massa è normalizzata sull’intervallo `[-k sigma,k sigma]`. In tre dimensioni il fattore `eta` vale 1.0081 e in LAPMOD è posto uguale a 1.

Kernel uniforme (CCA2)

Un'espressione generale per il calcolo delle concentrazioni in un punto `(x,y,z)` mediante Kernel è data da (Uliasz, 1994):

`C_(IST)(x,y,z) = sum_(i=1)^N m_i / (h_(x i) h_(yi) h_(zi)) [K(r_x,r_y,r_z)+K(r_x,r_y,r'_z)] `

dove:

• `N` è il numero di particelle computazionali che contribuiscono alla concentrazione nel punto di interesse
• `m_i` è la massa associata alla i-esima particella computazionale
• `h_(x i)`,`h_(yi)`,`h_(zi)` sono le larghezze di banda o bandwidths (m) della particella i-esima nelle tre direzioni

Le bandwiths determinano il grado di smoothing del campo di concentrazione; esse sono una funzione della posizione della particella e del tempo.

Le distanze adimensionali `r_x`, `r_y` e `r_z` sono definite come:

`r_x = (x_i - x) /h_x`

`r_y = (y_i - y) /h_y`

`r_z = (z_i - z) /h_z`

La funzione kernel K può essere qualsiasi a patto di soddisfare la condizione di normalizzazione:

`1 / (h_x h_y h_z) int_-oo^oo int_-oo^oo int_-oo^oo K dx dy dz = 1`

Per CCA2 sui recettori della griglia regolare, il profilo uniforme normalizzato è:

`K(r_x,r_y,r_z) = 1/8 I_x I_y I_z`

dove `I_i = 1` se `r_i^2 < 1`, altrimenti `I_i = 0`. Le larghezze di banda sono definite in maniera proporzionale alle griglie orizzontali di output (`Deltax, Deltay`) e all'altezza dello strato di rimescolamento (Hmix):

`h_x = alpha_x Deltax`

`h_y = alpha_y Deltay`

`h_z = alpha_z H_(mix)`

Un valore tipico per tre i fattori moltiplicativi `alpha_i` è 0.5 (Uliasz, 1994). Tali valori possono essere forniti in input in LAPMOD, quindi l'utente può valutare l'effetto della scelta dei fattori moltiplicativi sui valori di concentrazione calcolati. Valori elevati dei fattori di moltiplicazione conferiscono un maggiore smoothing al campo di concentrazione istantanea calcolato. Inoltre, maggiore è il valore dei fattori di moltiplicazione, maggiore è il numero di celle di output in cui una specifica particella computazionale contribuisce alla concentrazione. Piccoli valori dei fattori di moltiplicazione, tali da rendere le bandwiths molto minori delle ampiezze delle celle della griglia di output, rendono il metodo del kernel uniforme quasi simile al metodo del conteggio.

L’implementazione CCA2 per recettori discreti non include il corrispondente taglio verticale e il controllo della riflessione; non equivale quindi al kernel uniforme normalizzato su griglia. Questa limitazione permane dopo l’aggiornamento AGELIM. Vedere le limitazioni note.

Kernel parabolico (CCA3)

La funzione K è definita da Uliasz (1994):

`K(r_x,r_y,r_z) = 15 / (8 pi) (1 - r^2)I`

Le bande nominali sono calcolate indipendentemente dai coefficienti di input:

HX = min(HXA + HXB*sqrt(t) + HXC*t, BWMAX(1))
HY = min(HYA + HYB*sqrt(t) + HYC*t, BWMAX(2))
HZ = min(HZA + HZB*sqrt(t),         BWMAX(3))

t è l’età della particella in secondi e le bande sono in metri. HX e HY non devono necessariamente coincidere. Sono semiassi del supporto: in spazio libero la deviazione standard del kernel parabolico lungo ciascun asse è h/sqrt(7).

Rilascio graduale delle bande orizzontali

Per CCA2/3 con AGELIM > 0, dopo il calcolo delle bande nominali si applica il raccordo cubico AGELIM. Il limite iniziale è min(h, abs(OUTGRID)) e viene rilasciato progressivamente per 0 ≤ t < AGELIM. A partire da AGELIM si usano le bande nominali; HZ non cambia. Si mantiene il calcolo tramite kernel anche per particelle giovani, senza assegnarne l’intera massa a una sola cella.

code/omp_parallel/outputk.f90 calcola le bande effettive prima della normalizzazione e della scansione dei recettori. La regola conserva la dipendenza dalla dimensione della griglia e rende graduale il rilascio delle bande; non elimina la discontinuità al bordo del supporto di CCA2.

Consultare i dettagli di CCA1, CCA2 e CCA3 per profili, bande e limitazioni.

La concentrazione totale è la somma dei contributi delle particelle, eventualmente nulla:

`C(x,y,z) = sum_(i=1)^N C_i(x,y,z)`

N è il numero delle particelle nel dominio di simulazione.

Il programma principale invia CCA1 a OUTPUT e CCA2/3 a OUTPUTK. I valori sono valutati ai centri delle celle e ai recettori discreti. Entrambe le routine OpenMP usano accumuli per thread seguiti da una riduzione.

Il dominio di output è definito dalle coordinate degli angoli LL (inferiore sinistro) e UR (superiore destro). È possibile definire recettori discreti esterni alla griglia regolare. L’utente può specificare il numero di campionamenti tra due output consecutivi: la concentrazione viene valutata a ogni campionamento e l’output medio è ottenuto mediando i campionamenti.

Concentrazioni di picco (sostanze odorigene e infiammabili/esplosive)

L’impatto odorigeno in atmosfera viene solitamente valutato sulla base della concentrazione di odore espressa in unità odorimetriche al metro cubo (OUE/m3), che rappresentano il numero di diluizioni necessarie affinché il 50% degli esaminatori non avverta più l’odore del campione analizzato (norma UNI EN 13725:2004).

Quando LAPMOD è utilizzato per studi di impatto olfattivo, la concentrazione oraria di picco, che rappresenta la media su periodi brevi, ad esempio su 5 minuti o su 1 minuto, fino alla durata di un respiro (5 secondi), può essere calcolata in maniera classica applicando il peak-to-mean ratio (PMR) alle concentrazioni medie orarie in fase di postprocessamento mediante LAPOST. In questo caso il PMR è un fattore di moltiplicazione costante (pari a 2.3 secondo le Linee Guida della Regione Lombardia, 2012), che rende conto della fluttuazione dei valori della concentrazione oraria sul breve periodo (tipicamente 5 minuti) entro cui si manifesta il disturbo olfattivo.

Per condurre un calcolo secondo questo approccio semplificato è necessario effettuare una simulazione della dispersione di una sostanza qualsiasi e postprocessare i risultati mediante LAPOST, indicando in input RMULT = 2.3.

Il fattore 2.3 di cui sopra è una approssimazione adottata dalle Linee Guida per semplificare i calcoli e garantire un risultato conservativo che tuttavia può essere eccessivamnte penalizzante.

Mediante LAPMOD è anche possibile calcolare correttamente il peak to mean ratio secondo la formulazione originale di Smith (1973) e applicando la funzione di attenuazione esponenziale di Mylne (1992). Tale approccio è ad esempio quello implementato dal modello austriaco di dispersione degli odori AODM (Schauberger et al, 2000) e qui di seguito descritto.

Smith (1973) propone la seguente relazione per il peak-to-mean ratio:

`Psi_0 = C_p / C_m = ( t_m / t_p)^u`

dove `C_m` è la concentrazione mediata su un tempo lungo `t_m`, e `C_p` è la concentrazione di picco su un tempo breve `t_p`. Smith indica i seguenti valori di `u` basati sulle condizioni di stabilità atmosferica, individuate tramite le classi di stabilità di Paquill-Gifford:

Classe di stabilità	u
A	0.65 (*)
B	0.65
C	0.52
D	0.35
E	0.
F	0.

(*) In LAPMOD, nel caso di stabilità in classe A viene assunto il valore di B, non essendo presente tra i valori indicati da Smith (1973).

La scelta di `t_p = 5 s`, brevi quanto un respiro, si basa sulle misure di Mylne (1990). Questi valori sono validi solo in prossimità della sorgente, infatti grazie al rimescolamento turbolento il peak-to-mean ratio si riduce col trascorrere del tempo dal rilascio in atmosfera. Per descrivere ciò, Mylne e Mason (1991) hanno derivato dall'analisi delle fluttuazioni delle concentrazioni in un pennacchio la seguente espressione per il peak-to-mean ratio che include una attenuazione esponenziale dipendente dal rapporto tra tempo di volo e tempo Lagrangiano:

`Psi = 1 + (Psi_0 -1) \ exp (-0.7317 T / T_L)`

In LAPMOD la stabilità atmosferica viene individuata al momento del rilascio della particella e viene utilizzata per assegnare alla particella il valore `Psi_0`. Il tempo Lagrangiano viene calcolato come media pesata sugli spostamenti della particella nel tempo di volo `T = sum_i Delta t_i`, sulla base dei valori incontrati dalla particella lungo la propria traiettoria, e tale valore viene utilizzato per il calcolo del rapporto `T / T_L`:

`T/T_L = T^2 / (sum_i (T_{Li} Delta t_i))`

Il tempo Lagrangiano può essere interpretato come un indicatore della correlazione tra due tempi successivi del moto turbolento di una particella. Poiché il tempo Lagrangiano ha un valore diverso per le tre componenti spaziali, si utilizza cautelativamente il valore maggiore.

Il periodo di media `t_m` corrisponde all'intervallo di output delle concentrazioni di LAPMOD che l'utente specifica nel file di input e il tempo breve `t_p` è un parametro indicato anch'esso esplicitamente nel file di input. Si noti che per applicare questo metodo di calcolo è necessario che la concentrazione in output sia un valore medio ottenuto con più campionamenti. Deve essere quindi impostato in input NSAM > 1.

In maniera alternativa, la concentrazione oraria di picco può essere anche calcolata memorizzando la massima concentrazione stimata in ognuno dei periodi di campionamento che porteranno al calcolo della media oraria. Ad esempio, se la media oraria è calcolata utilizzando 12 campionamenti, verrà memorizzato il massimo valore medio di 5 minuti all'interno di ogni ora.

Il postprocessore LAPOST calcola inoltre i parametri FIDOL (Frequenza, Intensità, Durata, Offensività e Luogo), ad eccezione dell'offensività, che dipende dalla miscela odorigena rilasciata e può avere caratteristiche soggettive. Per i dettagli di calcolo si rimanda al manuale di LAPOST.

Riferimenti

• Bellasio R. and Bianconi R. (2026) Implementation of Building Downwash in the LAPMOD Atmospheric Dispersion Model. Chemical Engineering Transactions, Vol. 127, 241-246. doi:10.3303/CET26127041
• Baklanov A. and Sorensen J.H. (2001) Parameterisation of radionuclide deposition in atmospheric Long-Range transport modeling. Phys. Chem. Earth (B), Vol. 26, N. 10, pp. 787-799.
• Franzese P., Luhar A.K., Borgas M.S. (1999) An efficient Lagrangian stochastic model of vertical dispersion in the convective boundary layer. Atmospheric Environment, 2337-2345.
• Du S. (1997) Universality of the Lagrangian velocity structure function constant (` C_0`) across different kinds of turbulence. Boundary-layer Meteorology, 83, 207-219.
• Hanna S.R., Briggs G.A. and Hosker R.P. (1982) Handbook on atmospheric diffusion. Technical Information Center, U.S. Department of Energy, pp. 102.
• Hibberd M.F. and B.L. Sawford (1994) A saline laboratory model of the planetary boundary layer. Boundary-layer meteorology, 67, 229-250.
• Hurley P. and Physick W. (1993) A skewed homogeneous Lagrangian particle model for convective conditions. Atmospheric Environment 27A, 619-624.
• Janicke U. and Janicke L. (2001) A three-dimensional plume rise model for dry and wet plumes. Atmospheric Environment, Vol. 35, N. 5, 887-890.
• Knoll G.F. (1979) Radiation, detection and measurement. John Wiley and Sons, Inc. pp. 816.
• Luhar A.K. and Britter R. (1989) A random walk model for dispersion in inhomogeneous turbulence in a convective boundary layer. Atmospheric Environment, 23, 1911-1924.
• Manins P.C. (1979) Partial penetration of an elevated inversion layer by chimney plumes. Atmospheric Environment, 13, 733-741.
• Maryon R.H. and Ryall D.B. (1996) Developments to the UK nuclear accident response model (NAME). Department of Environment, UK Meteorological Office, DoE Report # DOE/RAS/96.011.
• Mylne K.R. (1990) Concentration fluctuation measurements of a tracer plume up to 1 km range in the atmosphere. Ninth Symposium on Turbulence and Diffusion. Roskilde, 168-171.
• Mylne K.R. (1992) Concentration fluctuation measurements in a plume dispersing in a stable surface layer. Boundary-layer Meteorology, 60, 15-48.
• Mylne K.R. (1990) Concentration fluctuation measurements of a tracer plume up to 1 km range in the atmosphere. Ninth Symposium on Turbulence and Diffusion. Roskilde, 168-171.
• Orville H.D., Hirsch J.H. and May L.E. (1980) Application of a cloud model to cooling tower plumes and clouds. Journal of Applied Meteorology, Vol. 19, 1260-1272.
• Presotto L., R.Bellasio and R.Bianconi (2005) Assessment of the visibility impact of a plume emitted by a desulphuration plant. Atmospheric Environment, 39(4), 719-737.
• Rezacova D. and Sokol S. (2000) On the influence of cooling towers on weather and climate. Czech Hydrometeorological Institute, 77 pages.
• Scire J.S., Robe F.R., Fernau M.E. and Yamartino R.J. (1999a) A user's guide for the CALMET meteorological model (Version 5.0). Earth Tech Inc., September 1999.
• Scire J.S., Strimaitis D.G., and Yamartino R.J. (1999b) A user's guide for the CALPUFF dispersion model (Version 5.0). Earth Tech Inc., May 1999.
• Seinfeld J.H. and S.N. Pandis (1998) Atmospheric chemistry and physics. From air pollution to climate change, John Wiley & Sons.
• Shauberger G., M. Piringer and E. Pets (2000) Diurnal and annual variation of the sensation distance of odour emitted by livestock buildings calculated by the Austrian odour dispersion model (AODM). Atmospheric Environment, 34, 4839-4851.
• Smith M.E. (1973) Recommended Guide for the Prediction of the Dispersion of Airborne Effluents. ASME, New York.
• Thomson D.J. (1987) Criteria for the selection of stochastic models of particle trajectory in turbulence flows. J.Fluid Mech., 180, 529-556.
• Uliasz, M. (1994) Lagrangian particle dispersion modeling in mesoscale applications, in: Environmental Modeling, Vol. II, edited by: Zannetti, P., Computational Mechanics Publications, Southampton, UK.
• Verver G and F.A.De Leeuw (1992) An operational puff dispersion model. Atmospheric Environment, 26,17,3179-3193.
• Webster H.N. and Thomson D.J (2002) Validation of a Lagrangian model plume rise scheme using the Kincaid data set. Atmospheric Environment, Vol. 36, N. 32, pp. 5031-5042.
• Weil J. C. (1990) A diagnosis of the asymmetry in top down and bottom up diffusion using a Lagrangian stochastic model. Journal of Atmospheric Science, 47, 501-515.
• Yamada T. and Bunker S. (1988) Development of a nested grid, second moment turbulence closure model and application to the 1982 ASCOT Brush Creek data simulation. Journal of Applied Meteorology, 27, 562-578.
• Zhang L., Gong S., Padro J. and Barrie L. (2001) A size segregated particle dry deposition scheme for an atmospheric aerosol module. Atmospheric Environment, 35, 549-560.
• Petersen R.L. and Guerra S.A. (2018) PRIME2: Development and evaluation of improved building downwash algorithms for rectangular and streamlined structures. Journal of Wind Engineering and Industrial Aerodynamics, 173, 67-78. doi:10.1016/j.jweia.2017.11.027
• Petersen R.L., Paumier J.O. and Guerra S.A. (2022) Development, evaluation, and implementation of building downwash and plume rise enhancements in AERMOD. Journal of the Air & Waste Management Association, 72(12), 1423-1441. doi:10.1080/10962247.2022.2120563
• Schulman L.L., Strimaitis D.G. and Scire J.S. (2000) Development and Evaluation of the PRIME Plume Rise and Building Downwash Model. Journal of the Air & Waste Management Association, 50(3), 378-390. doi:10.1080/10473289.2000.10464017
• Wilson D.J. (1979) Flow patterns over flat-roofed buildings and application to exhaust stack design. ASHRAE Transactions, 85(2), 284-295.