STIMA DELLA PROFONDITà DELLA MOHO NEI BALCANI OCCIDENTALI DA OSSERVAzIONI DI GRAVITà DEL SATELLITE GOCE D. SampietroGReD s.r.l., Como

Introduzione. L’area dei Balcani occidentali, ovvero la zona che si estende tra la Bulgaria e il Mar Adriatico, è una delle regioni europee più complesse e attive dal punto di vista tettonico. L’area si trova in corrispondenza della collisione tra la placca africana e quella eurasiatica ed è caratterizzata dalla presenza della cintura orogenetica alpino-himalayana e dall’apertura del bacino pannonico. Infatti, la spinta della placca Adriatica nella litosfera europea, ha causato la formazione di importanti catene montuose come le Alpi, le Dinaridi e le Albanidi. Secondo studi più recenti, basati principalmente sull’analisi delle velocità da reti GNSS permanenti, e da metodi di sismica passiva la placca Adriatica può essere suddivisa in una serie di due o anche tre unità più piccole in movimento verso nord con velocità dell’ordine di 35 mm/anno (Herak et al., 2005; Ivančić et al., 2006). Al contrario la parte meridionale dei Carpazi e la parte orientale della penisola balcanica mostrano un movimento orientato verso sud di circa 3 mm/anno. Oltre che per l’estrema complessità geologica i Balcani occidentali rappresenta anche una zona interessante anche da un punto di vista “storico”: infatti i primi studi sulla litosfera e in particolare sulla determinazione della Moho sono stati condotti in questa regione da Mohorovičić che per primo, studiando il terremoto avvenuto nella valle di Kupa nel 1909 (Mohorovičić, 1992), individuò la presenza di una discontinuità tra crosta e mantello terrestre, la cosiddetta discontinuità di Mohorovičić o Moho. Nel suo lavoro Mohorovičić, sulla base delle registrazioni di una serie di terremoti, osservava la presenza di due distinte coppie di onde P e S una delle quali generata a una discontinuità strutturale sotto la superficie della Terra. Nel suo articolo Mohorovičić stimò la profondità di questa discontinuità nella zona croata in circa 54 km. Da allora sono stati condotti molti altri studi per comprendere meglio la struttura litosferica sotto i Balcani occidentali: a partire dai lavori di Dragaševií e Andric (1968) che utilizzarono due profili di sismica profonda lungo le Alpi Dinariche, al lavoro di Aljinović et al. (1984) che suggerisce, sulla base di altri tre profili sismici che vanno dalla costa adriatica al continente in direzione sud est-nord ovest, uno spessore della crosta di 45 km sotto le Dinaridi rapidamente decrescente fino ad arrivare a soli 20 km nel bacino pannonico e in corrispondenza del mar Adriatico.

In aggiunta è importante ricordare qui l’analisi della velocità di propagazione delle onde sismiche nella regione circum adriatica di Herak e Herak (1995) che ha mostrato uno spessore crostale media di 40 km nelle Dinaridi crescente verso sud-est e fino a raggiungere 55 km nella loro parte più meridionale, il lavoro di Van der Meijde et al. (2003) e di Stipčević et al. (2011) che hanno analizzato i tempi d’arrivo delle onde sismiche in corrispondenza di stazioni sismologiche per stimare modelli semplici della struttura della crosta terrestre.

In particolare in Van der Meijde et al. (2003) due stazioni lungo la costa croata sono stati indagate stimando uno spessore crostale di 47 e 41 km, con una incertezza dell’ordine di 1.6 km mentre in Stipčević et al. (2011), una serie di 8 stazioni è stata utilizzata per trovare uno spessore medio della crosta sotto le stazioni dell’Adriatico settentrionale in linea con i risultati pubblicati dai recenti esperimenti DSS [profilo Alp07: Šumanovac et al. (2009)], mentre nelle Dinaridi centrali e nel nel sud della Croazia hanno osservato una Moho significativamente più profonda delle stime precedenti.

È importante notare che tutti questi studi derivanti da osservazioni sismiche soffrono ancora di mancanza di osservazioni nella regione delle Dinaridi e nelle aree circostanti, che si traduce inevitabilmente in elevate incertezze nelle stime della profondità Moho. Infatti, mentre la struttura crostale del bacino pannonico è ragionevolmente ben conosciuta grazie all’esplorazione delle compagnie petrolifere, si veda ad esempio Dolton (2006) la disponibilità di osservazioni in Croazia e nella vicina Bosnia-Erzegovina è abbastanza limitata.

1��

GNGTS 2014 sessione 3.2

001-258 volume 3 196 24-10-2014 16:46:52

GNGTS 2014 sessione 3.2

1��

Nel presente lavoro i dati di gravità, e in particolare la quinta release del modello di gravità globale ottenuto da osservazioni GOCE applicando il metodo Time-Wise (Pail et al., 2010) sono stati utilizzati per dedurre informazioni sulla struttura crostale e sulla profondità della Moho nei Balcani occidentali. La procedura può essere divisa in due fasi principali: la prima consiste nel riconoscere e isolare le diverse province geologiche nell’area di studio sfruttando le informazioni provenienti dal modello globale del campo gravitazionale stesso, mentre la secondo consiste nell’invertire il campo gravitazionale per stimare la profondità della Moho e alcune informazioni sulla densità della crosta.

Nel secondo capitolo l’algoritmo e i risultati della classificazione in province geologiche saranno descritti e discussi, mentre nel terzo capitolo verrà presentato l’algoritmo di inversione. Infine nell’ultimo capitolo verranno esposti e analizzati i principali risultati numerici.

Definizione delle province geologiche. Per poter stimare la profondità della Moho da osservazioni del campo gravitazionale è necessario conoscere le principale variazioni di densità all’interno della crosta e nei primi strati del mantello. Infatti una volta che tali variazioni sono note è possibile calcolarne l’effetto in termini di campo gravitazionale e quindi isolare e rimuovere dai dati il segnale dovuto alla discontinuità di Mohorovičić. Per questo è necessario modellizzare la geometria e la densità dei sedimenti, della crosta cristallina, di eventuali mari o ghiacci presenti e del mantello superiore. Per quanto riguarda la crosta cristallina, una possibilità per modellizzare almeno le principali variazione di densità laterali, studiata e applicata recentemente a livello globale in Reguzzoni e Sampietro (2014), consiste nel suddividere la crosta in regioni geologicamente omogenee, ognuna delle quali classificate come una di otto tipi di crosta (i.e. scudi, piattaforma continentale, bacini sedimentari, crosta in estensione, zone orogenetiche, crosta oceanica, ridge oceanici e province ignee). Per ogni classi di crosta viene quindi definita una funzione empirica [basata sul lavoro di Christensen and Mooney (1995)] che descriva variazioni di densità rispetto a variazione di profondità. È importante notare che tali funzioni sono richieste oltre che per la riduzione dei dati anche per definire il contrasto di densità tra crosta e mantello terrestre necessario per l’inversione del segnale gravitazionale residuo.

In questo contesto la definizione dei limiti geografici delle province geologiche è un tema cruciale per ridurre correttamente i dati e invertire il segnale residuo. Nel presente lavoro è proposto un approccio bayesiano per classificare a partire da un funzionale del campo gravitazionale, la regione indagata in province geologiche.

In dettaglio consideriamo una griglia di anomalie di gravità δg ad una certa altitudine h costante, possibilmente vicina alle masse topografiche, ottenuta applicando un operatore di

Fig. 1 – Mappa delle province geologiche a-priori digitalizzata a 1°x1° (a) e risultato dell’algoritmo di classificazione proposto (b).

001-258 volume 3 197 24-10-2014 16:46:56

sintesi armonica (Moritz, 1980) a un modello di campo gravitazionale globale e supponiamo di conoscere un modello approssimato, ad esempio da una mappa a-priori delle province geologiche (Exxon, 1995), di province geologiche Gi per ogni nodo i della griglia. In generale, Gi può assumere n valori dove n è il numero delle province geologiche nella zona di interesse. Si veda ad esempio la mappa delle province geologiche dei Balcani occidentali con una risoluzione di 1°x1° ottenuta digitalizzando la mappa delle province geologiche sviluppata dall’U.S. Geological Survey riportata in Fig. 1. Dal momento che solo quattro province geologiche sono presenti, n = 4, Gi può assumere solo un valore nel set delle province disponibili, i.e. : G1, G2, G3, G4.

Per ogni pixel i inoltre viene selezionato un intorno Δi, ad esempio formato dagli 8 pixel adiacenti a quello considerato. La probabilità a priori che Gi assuma un certo valore Gk con k = 1,2,…,n, cioè P(Gi = Gk) è calcolata per mezzo di una matrice peso W:

dove (1)

La matrice dei pesi W è definita in base alla distanza tra ciascun pixel in j = {i,Δi} e il pixel i. Ad esempio, supponendo che in un intorno di un certo pixel i il modello a priori abbia solo due province geologiche con la geometria e la matrice dei pesi W di Fig. 2 le due probabilità a priori sono e .

Fig. 2 – Esempio di matrice dei pesi. Il pixel i è posto nel centro della matrice.

Per calcolare la likelihood, i.e. P(δgi|Gi = Gk) supponiamo che la crosta sia in perfetto equilibrio isostatico (secondo il modello di Airy) è che l’effetto gravitazionale ad ogni nodo della griglia sia in prima approssimazione dovuto ad una slab di Bouguer con spessore pari allo spessore della crosta. Considerando queste approssimazioni è possibile scrivere una relazione lineare tra δgi e la densità della crosta ρGii

. A questo punto, partendo dal modello a-priori di province geologiche è possibile stimare ai minimi quadrati ρ̂Gii

per ogni provincia geologica e la relativa varianza σ̂ 2ρGii. Infine per ogni pixel la likelyhood è calcolata come la probabilità che la densità di un certo pixel i appartenga a una distribuzione normale con media ρ̂Gii

e varianza σ̂ 2ρGii. La probabilità a posteriori è calcolata applicando il noto teorema di Bayes come il prodotto tra la probabilità a-priori e la likelihood. La provincia geologica del singolo pixel i è quindi scelta in modo da massimizzare la probabilità a posteriori.

A questo punto due osservazioni sono necessarie: in primo luogo è importante sottolineare che le approssimazioni introdotta dall’ipotesi di crosta perfettamente isostatica e dal calcolo del segnale gravitazionale utilizzando la semplice plate di Bouguer non consentono una corretta determinazione della densità crostale, tuttavia questa prima fase ha come scopo finale solamente quello di individuare regioni omogenee (dal punto di vista della struttura della crosta) all’interno dell’area di studio. In secondo luogo va anche notato che l’uso della plate di Bouguer permette di considerare Gi come una realizzazione di un campo di Markov con correlazione significativa solo con i vicini più prossimi permettendo l’applicazione di un metodo come il Gibbs sampler per massimizzare la probabilità a posteriori (Smith e Roberts, 1993; Sansò et al., 2011) e quindi facilitando notevolmente la classificazione delle zone omogenee.

Il metodo è stato applicato per migliorare la modellazione dei principali confini delle province geologiche nei Balcani occidentali. In particolare è stato utilizzato come punto di partenza dell’algoritmo di classificazione il modello a priori mostrato in Fig. 1, considerando come osservazioni una griglia di gravity disturbances sintetizzate dal modello globale GO_CONS_GCF_2_TIM_R5 (Pail et al. 2010) ad una quota di 3000 m con risoluzione di 1’ e

1��

GNGTS 2014 sessione 3.2

001-258 volume 3 198 24-10-2014 16:46:57

GNGTS 2014 sessione 3.2

1��

considerando una matrice peso W (come in Fig. 2) definita da α = 0.3 e β = γ = (1 – α)/8.Va sottolineato come questa configurazione della matrice dei pesi obblighi l’algoritmo a

scegliere la provincia geologica per il pixel i solo tra le province geologiche adiacenti andando quindi a modificare solamente i confini tra due (o più) province.

Il risultato finale è mostrato in Fig. 1 dove si può vedere come il metodo sia in grado di modificare correttamente i confini delle province geologiche ottenendo una mappa molto più dettagliata rispetto a quella di partenza e sostanzialmente guidata solamente dalle osservazioni del campo gravitazionale.

Stima della Moho. L’algoritmo di inversione si basa sulla soluzione locale sviluppata all’interno del progetto GEMMA (GOCE Exploitation for Moho Modeling and Applications) finanziato dal programma dell’Agenzia Spaziale Europea Support To Science Element (STSE). Nel seguito sono riportati solamente i concetti principali, rimandando il lettore a consultare Sampietro (2011), Reguzzoni e Sampietro (2012) e Sampietro et al. (2014) per dettagli.

Prima di descrivere l’algoritmo di inversione è importante sottolineare il sistema di riferimento utilizzato nel presente studio. Essendo l’area di interesse abbastanza ridotta la soluzione si basa su un’approssimazione planare del problema. Infatti è stato dimostrato che per regioni con estensione inferiore a 10°x10° la differenza tra l’approssimazione planare e l’approssimazione sferica è trascurabile essendo inferiore a 0.5 km in termine di profondità della Moho (Sampietro, 2011). La prima operazione consiste quindi nel mappare il sistema di riferimento globale geodetico (i.e. latitudine, longitudine e quota ellissoidica) in un sistema di coordinate cartesiane locali. Questo mapping è definito dalle seguenti equazioni:

(2)

dove ϕ, λ e h sono la latitudine, longitudine e altezza ellissoidica, rispettivamente, di un certo nodo della griglia delle osservazioni, ϕ– e λ– e sono la latitudine e la longitudine del centro della regione considerata, R è il raggio della sfera locale e infine x, y, e z sono le coordinate mappate. Si noti che questa operazione è solo un cambiamento di coordinate e che il nuovo sistema di riferimento è solo approssimativamente quello definito da una terna locale con origine nel centro della griglia e tangente all’ellissoide (Sansò, 2006).

Considerando il sistema di riferimento di cui sopra, è possibile applicare l’algoritmo di inversione basato su un processo di deconvoluzione di Wiener nel dominio delle frequenze che filtra gli errori delle osservazioni e rende numericamente efficiente la risoluzione del problema [per dettagli si veda Reguzzoni e Sampietro (2012)]. L’unicità della soluzione è garantita dall’approssimazione della struttura crosta-mantello con un semplice modello a due strati (Sampietro e Sansò, 2012), quest’approssimazione è resa valida rimuovendo, a monte dell’inversione del campo gravitazionale, gli effetti dovuti alle principali variazioni di densità (e.g. topografia, batimetria, sedimenti e mantello superiore). La soluzione è stata quindi migliorata adattando la strategia globale proposta in Reguzzoni e Sampietro (2014) all’inversione locale. In sintesi tale miglioramento permette, conoscendo informazioni derivanti ad esempio da profili sismici, di stimare la profondità media della Moho, un fattore di scala per la funzione profondità-densità di ogni provincia geologica e di considerare nella soluzione eventuali variazioni di profondità della densità della crostale.

Per quanto riguarda la riduzione del segnale gravitazionale è stato utilizzato il modello ETOPO1 (Amante e Eakins, 2009) per la topografia e la batimetria, un modello di sedimenti a 1°x1° (Laske e Master, 1997) e il modello GyPSuM per il mantello superiore (Simmons et al., 2010). Come informazione sismica per la stima della profondità media della Moho è stato utilizzato il modello CRUST1.0 (Laske et al., 2013 ). I risultati dell’inversione in termini di profondità della Moho sono presentati in Fig.3

001-258 volume 3 199 24-10-2014 16:46:57

Conclusioni. La Moho stimata sembra presentare tutte le caratteristiche note in letteratura quali l’elevato spessore in corrispondenza delle maggiori regioni orogenetiche e la riduzione dello spessore della crosta in corrispondenza della placca adriatica e dei bacini sedimentari. Lo spessore della crosta in corrispondenza della Pannonia, compreso tra 20 e 30 km è consistente con diversi modelli di crosta (e.g. Grad e Tira, 2009). Questi risultati mostrano in primo luogo la capacità dell’osservazioni della missione GOCE di classificare regioni omogenee di crosta terrestre anche in presenza di geologie complesse. Inoltre dimostrano la bontà dell’algoritmo di inversione utilizzato e l’elevata accuratezza e risoluzione che si può ottenere sfruttando

localmente le osservazioni globali del campo gravitazionale.

BibliografiaAmante, C., & Eakins, B. W. (2009). ETOPO1 1 arc-minute global relief model: procedures, data sources and

analysis (p. 19). US Department of Commerce, National Oceanic and Atmospheric Administration, National Environmental Satellite, Data, and Information Service, National Geophysical Data Center, Marine Geology and Geophysics Division.

Aljinović, B., Blašković, I., Cvijanović, D., Prelogović, E., Skoko, D., & Brdarević, N. (1984). Correlation of geophysical, geological and seismological data in the coastal part of Yugoslavia. Boll. Ocean. Teor. Appl, 2, 77-90.

Christensen, N. I., & Mooney, W. D. (1995). Seismic velocity structure and composition of the continental crust: A global view. Journal of Geophysical Research: Solid Earth (1978–2012), 100(B6), 9761-9788.

Dolton, G. L. (2006). Pannonian Basin Province, Central Europe (Province 4808)--petroleum Geology, Total Petroleum Systems, and Petroleum Resource Assessment. US Department of the Interior, US Geological Survey.

Dragaševií, T., & Andric, B. (1968). Deep Seismic Sounding of the Earth’s Crust in the Area of the Dinarides and the Adriatic Sea. Geophysical Prospecting,16(1), 54-76.

Exxon, 1995. Tectonic Map of the World, 18 sheets, scale 1: 10,000,000. 594 Technical Report. Exxon, Houston , Texas.

Grad, M., & Tiira, T. (2009). The Moho depth map of the European Plate. Geophysical Journal International, 176(1), 279-292.

Herak, D., & Herak, M. (1995). Body-wave velocities in the circum-Adriatic region. Tectonophysics, 241(1), 121-141.

Herak, D., Herak, M., Prelogović, E., Markušić, S., & Markulin, Ž. (2005). Jabuka island (central Adriatic Sea) earthquakes of 2003. Tectonophysics, 398(3), 167-180.

Ivančić, I., Herak, D., Markušić, S., Sović, I., & Herak, M. (2006). Seismicity of Croatia in the period 2002-2005. Geofizika, 23(2), 87-103.

Laske G. & Masters G. (1997). A Global Digital Map of Sediment Thickness, EOS Trans. AGU, 78, F483.Laske, G., Masters, G., Ma, Z., & Pasyanos, M. (2013). Update on CRUST1.0 - A 1-degree Global Model of Earth’s

Crust, Geophys. Res. Abstracts, 15, Abstract EGU2013-2658.Mohorovičić, A. (1992). Earthquake of 8 October 1909. Geofizika, 9(1), 3-55.Moritz, H. (1980). Advanced physical geodesy. Karlsruhe: Wichmann; Tunbridge, Eng.: Abacus Press, 1980. 1.Pail, R., Goiginger, H., Mayrhofer, R., Schuh, W. D., Brockmann, J. M., Krasbutter, I., & Fecher, T. (2010). GOCE

gravity field model derived from orbit and gradiometry data applying the time-wise method. In Proceedings of the ESA Living Planet Symposium, ESA Publication SP-686, ESA/ESTEC, ISBN (Online) (pp. 978-92).

Reguzzoni, M., & Sampietro, D. (2012). Moho estimation using GOCE data: a numerical simulation. In Geodesy for Planet Earth (pp. 205-214). Springer Berlin Heidelberg.

Reguzzoni, M., & Sampietro, D. (2014). GEMMA: An Earth crustal model based on GOCE satellite data. International Journal of Applied Earth Observation and Geoinformation.

Fig. 3 – Stima della profondità della Moho (in km) nell’area di studio.

200

GNGTS 2014 sessione 3.2

001-258 volume 3 200 24-10-2014 16:46:59

GNGTS 2014 sessione 3.2

201

Smith, A. F., & Roberts, G. O. (1993). Bayesian computation via the Gibbs sampler and related Markov chain Monte Carlo methods. Journal of the Royal Statistical Society. Series B (Methodological), 3-23.

Sampietro, D. (2011). GOCE exploitation for Moho modeling and applications. In Proc. of the 4th International GOCE User Workshop (Vol. 31).

Sampietro, D., Reguzzoni, M., & Braitenberg, C. (2014). The GOCE estimated Moho beneath the Tibetan Plateau andThe GOCE estimated Moho beneath the Tibetan Plateau and Himalaya. In Earth on the Edge: Science for a Sustainable Planet (pp. 391-397). Springer Berlin Heidelberg.Springer Berlin Heidelberg.

Sansò, F. (2006). Navigazione geodetica e rilevamento cinematico. Polipress.Polipress.Sampietro, D., & Sansó, F. (2012). Uniqueness theorems for inverse gravimetric problems. In VII Hotine-Marussi

Symposium on Mathematical Geodesy (pp. 111-115). Springer Berlin Heidelberg.Simmons, N. A., Forte, A. M., Boschi, L., & Grand, S. P. (2010). GyPSuM: A joint tomographic model of mantle

density and seismic wave speeds. Journal of Geophysical Research: Solid Earth, 115(B12).Stipčević, J., Tkalčić, H., Herak, M., Markušić, S., & Herak, D. (2011). Crustal and uppermost mantle structure

beneath the External Dinarides, Croatia, determined from teleseismic receiver functions. Geophysical journal international, 185(3), 1103-1119.

Šumanovac, F., Orešković, J., & Grad, M. (2009). Crustal structure at the contact of the Dinarides and Pannonian basin based on 2-D seismic and gravity interpretation of the Alp07 profile in the ALP 2002 experiment. Geophysical Journal International, 179(1), 615-633.

Van Der Meijde, M., Van Der Lee, S., & Giardini, D. (2003). Crustal structure beneath broad-band seismic stations in the Mediterranean region. Geophysical Journal International, 152(3), 729-739.

001-258 volume 3 201 24-10-2014 16:47:00