Laboratorio Metodi di Potenziale Prof Carla Braitenberg Tutor

Laboratorio Metodi di Potenziale Prof. Carla Braitenberg Tutor Dott. Tommaso Pivetta AA. 2018 -2019 Università di Trieste

Install GMT 5. 0 • You must install GMT on your PC for the laboratory. At the course windows system is used. If your system is different from Windows, please use the PCs of the University or install a Windows partition on your pc. • GMT is a versatile tool for mapping and modeling the potential field data. We will use it to display the gravity and magnetic fields. • Download link: http: //gmt. soest. hawaii. edu/projects/gmt/wiki/Download • The mapping tool is a collection of single routines which are called from the command line. They read input data and create an output plot file (default is postscript). To display the plot file you must launch a secondary program that reads the postscript file. • gmt-5. 4. 4 -darwin-x 86_64. dmg (installer for Mac) • gmt-5. 4. 4 -win 64 (installer for Windows) • gmt-5. 4. 4 -win 32 (installer for Windows)

GMT other useful links • Link al software: http: //gmt. soest. hawaii. edu/projects/gmt/wiki/Installing • Link alla documentazione: • http: //gmt. soest. hawaii. edu/doc/5. 4. 3/index. html • E’ consigliato installare un programma per la lettura dei file in formato Postscript: gv (ftp: //ftp. gnu. org/gnu/gv/) o ghostview • Anche un buon editor di testo, come Notepad++ • https: //qpdownload. com/notepad-plus/

The files offered for download

Double click on the win 32 or win 64. exe file

By default the GMT 5. 0 is installed under the directory programs

You need a postscript viewer. GMT recommends Post. Script previewer (e. g. , gv (or ghostview or plain ghostscript)) Once the ghostview has installed (for instance in C: Programmi), please copy the folder gs to the folder c: programs

The documentation and examples of GMT are stored in the folder doc • C: programsgmt 5sharedoc • The subfolder html contains a tutorial and the full instructions on how to use the routines • Double-click on GMT_Tutorial. html for a tutorial or GMT_Docs for a full documentation.

To run the GMT routines open a command line window. You can type cmd.

By typing the command gmt on the command line the below text is diplayed. It means that GMT was installed correctly on your pc.

It is necessary that your PC has the correct international settings of how to display numbers

You must set the decimal separator to dot, not comma. Here the default settings, which are not ok for GMT.

After you have corrected the settings you should have the following:

For launching the GMT commands, save the commands to a script. Save all scripts in the folder C: MEPO. Please create a working folder C: MEPOinput_data. To be executable the script must end with the extension myscript. bat. Double click on the bat-file will launch the script. The following script will create the map you see in the next slide. Install the program Notepad++ for editing the scripts. • gmt pscoast -R-90/-70/0/20 -JM 6 i -P -Ba -Gchocolate > GMT_tut_3. ps • c: Programsgsgs 9. 25bingswin 64 GMT_tut_3. ps

From the GMT tutorial • gmt pscoast -R-90/-70/0/20 -JM 6 i -P -Ba -Gchocolate > GMT_tut_3. ps • c: Programsgsgs 9. 25bingswin 64 GMT_tut_3. ps

Explanation of the pscoast command. • pscoast -R-5/10/-5/5 -JM 6 i -P -Ba -Gchocolate > Mepo_Es 1. ps • -R-5/10/-5/5: estremi in longitudine e tatitudine • -JM 6 i: proezione di mercatore. Mappa larga 6 inch • -P: portrait (orientazione della carta) • -Ba: annotazione degli assi • -Gchocolate: terre emerse in beige • > Mepo_Es 1. ps: output su file di tipo post-script

Mappa di Coste con coordinate • psbasemap -R-5/10/-5/5 -JM 6 i -P -Ba -B+t"Esercizio 2" -K > Mepo_Es 2. ps • Psbasemap: mappa di base con le coordinate • > apri file di uscita • >> continua a scrivere sullo stesso file • -K: arrivera’ altro materiale per il file • pscoast -R-5/10/-5/5 -JM 6 i -P –Ba -Gchocolate -O >> Mepo_Es 2. ps • -O: overlay sul file precedente

Ulteriori approfondimenti: • Leggere la documentazione delle singole istruzioni sul sito GMT • Creare ulteriori mappe • Creare mappa globale delle coste, utilizzando altre proiezioni geografiche. Si veda il GMT tutorial 1: • http: //gmt. soest. hawaii. edu/doc/5. 3. 2/GMT_Tutorial. html#sessionone

Indice degli argomenti da trattare durante il laboratorio 1 Installazione e primo utilizzo del pacchetto software GMT (General Mapping Tools). Installare anche database di topografia e coste e fiumi. Creazione di mappe con topografia e coste. 2 Creazione di mappe di: anomalia di gravita’ , Bouguer e topografia per zona a scelta dello studente. Recupero dati dal database ICGEM 3 Modellazione di un corpo tramite un prisma. Analisi delle proprieta’ del campo in funzione delle dimensioni e profondita’ del corpo. Vedi esempi sul testo Hinze, von Freese, Saad, pagina 183. 4 Modellazione di una area tramite tesseroidi 5 Calcolo di anomalia isostatica

Recuperare i dati dal sito ICGEM • Le schede successive illustrano come creare i dati di gravita’ dal sito ICGEM e rappresentarli graficamente in GMT. • http: //icgem. gfz-potsdam. de/home

Calculate the gravity anomaly at zero height for a static gravity field model in spherical harmonics. Use model EIGEN 6 C 4, choose classic gravity anomaly. For definition of the functionals the user guide of the ICGEM page.

You must download the grid. You will find the grid in the download folder. Please copy it to your working folder c: MEPOdati_input. The next slide shows the header of the file. You can open the file in Notepad++.

The header contains all the information abut the spherical harmonic model used to make the calculations, the parameters used. Check the number of lines use by the header- the number of lines must be known when we intend to read the file in GMT.

Convert column format to grid format of GMT and plot • xyz 2 grd. input_dataEIGEN 6 C 4_grav_Anom_cl. gdf G. input_dataGrav. Anom_Cl. grd -h 34 -I 0. 1/0. 1 R 10/30/20/50 -V • grdinfo. input_dataGrav. Anom_Cl. grd latlimit_north 50. 000000 latlimit_south 20. 000000 longlimit_west 10. 000000 longlimit_east 30. 000000

Plot the grid as contour lines on top of the basemap • psbasemap -R 10/30/20/50 -JM 5 i -P -Ba -B+t"Gravity Anomaly" -K > Mepo_Es 1. ps • grd 2 cpt. input_dataGrav. Anom_Cl. grd -Crainbow -L-180/180 -Z > tess. cpt • pscoast -R 10/30/20/50 -JM 5 i -I 1 -Na -Ggray -Dh -P -Ba -O -K >> Mepo_Es 1. ps • grdcontour. input_dataGrav. Anom_Cl. grd -R 10/30/20/50 -JM 5 i -C 30 -O >> Mepo_Es 1. ps • c: Programsgsgs 9. 25bingswin 64 -s. DEVICE=jpeg -d. Use. Crop. Box -r 300 s. Output. File=Mepo_Es 1. jpg Mepo_Es 1. ps

Plot the grid as contour lines and color plot. Here the coastline is drawn. • xyz 2 grd. input_dataEIGEN-6 C 4_grav_Anom_cl. gdf -G. input_dataGrav. Anom_Cl. grd -h 34 -I 0. 1/0. 1 R 10/30/20/50 -V • grdinfo. input_dataGrav. Anom_Cl. grd • psbasemap -R 10/30/20/50 -JM 5 i -P -Ba -B+t"Gravity Anomaly" -K > Mepo_Es 1 b. ps • • grd 2 cpt. input_dataGrav. Anom_Cl. grd -Crainbow -L-180/180 -Z > tess. cpt grdimage -R 10/30/20/50. input_dataGrav. Anom_Cl. grd -Ctess. cpt -JM 5 i -O -K >> Mepo_Es 1 b. ps pscoast -R 10/30/20/50 -JM 5 i -I 1 -Na -W -Dh -P -Ba -O -K >> Mepo_Es 1 b. ps grdcontour. input_dataGrav. Anom_Cl. grd -R 10/30/20/50 -JM 5 i -C 30 -O >> Mepo_Es 1 b. ps • c: Programsgsgs 9. 25bingswin 64 -s. DEVICE=jpeg -d. Use. Crop. Box -r 300 -s. Output. File=Mepo_Es 1 b. jpg Mepo_Es 1 b. ps

Resulting maps of previous scripts

The ICGEM calculation service allows to calculate also other functionals. Calculate the Bouguer field (simple Bouguer), the geoid, and the topography for the same area. Leave limits and gravity model the same. Here the screen for the simple Bouguer anomaly.

Plot the Bouguer as contour lines on top of the basemap. • xyz 2 grd. input_dataEIGEN-6 C 4_grav_BG. gdf -G. input_dataGrav. Anom_BG. grd -h 37 -I 0. 1/0. 1 R 10/30/20/50 -V • grdinfo. input_dataGrav. Anom_BG. grd • psbasemap -R 10/30/20/50 -JM 5 i -P -Ba -B+t"Gravity Bouguer" -K > Mepo_Es 1_BG. ps • grd 2 cpt. input_dataGrav. Anom_BG. grd -Crainbow -L-220/320 -Z > tess. cpt • grdimage -R 10/30/20/50. input_dataGrav. Anom_BG. grd=sf -Ctess. cpt -JM 5 i -O -K >> Mepo_Es 1_BG. ps • pscoast -R 10/30/20/50 -JM 5 i -I 1 -Na -W -Dh -P -Ba -O -K >> Mepo_Es 1_BG. ps • grdcontour. input_dataGrav. Anom_Cl. grd=sf -R 10/30/20/50 -JM 5 i -C 30 -O >> Mepo_Es 1_BG. ps • c: Programsgsgs 9. 25bingswin 64 -s. DEVICE=jpeg -d. Use. Crop. Box -r 300 s. Output. File=Mepo_Es 1_BG. jpg Mepo_Es 1_BG. ps Attention, number of header lines is different!

create and plot the Geoid map • : : ------Geoid • xyz 2 grd. input_dataEIGEN-6 C 4_Geoid. gdf -G. input_dataGeoid. grd -h 36 -I 0. 1/0. 1 R 10/30/20/50 -V • grdinfo. input_dataGeoid. grd • psbasemap -R 10/30/20/50 -JM 5 i -P -Ba -B+t"Geoid" -K > Mepo_Es 1_GE. ps • grd 2 cpt. input_dataGeoid. grd -Crainbow -L-150/150 -Z > tess. cpt • grdimage -R 10/30/20/50. input_dataGeoid. grd -Ctess. cpt -JM 5 i -O -K >> Mepo_Es 1_GE. ps • pscoast -R 10/30/20/50 -JM 5 i -I 1 -Na -W -Dh -P -Ba -O -K >> Mepo_Es 1_GE. ps • grdcontour. input_dataGeoid. grd -R 10/30/20/50 -JM 5 i -C 30 -O >> Mepo_Es 1_GE. ps • c: Programsgsgs 9. 25bingswin 64 -s. DEVICE=jpeg -d. Use. Crop. Box -r 300 s. Output. File=Mepo_Es 1_GE. jpg Mepo_Es 1_GE. ps

Create topography in ICGEM and plot it • : : ------Etopo 1 • xyz 2 grd. input_dataetopo 1_selection. gdf -G. input_dataEtopo 1_select. grd -h 36 -I 0. 1/0. 1 R 10/30/20/50 -V • grdinfo. input_dataEtopo 1_select. grd • psbasemap -R 10/30/20/50 -JM 5 i -P -Ba -B+t"Etopo 1" -K > Mepo_Es 1_Topo. ps • grd 2 cpt. input_dataEtopo 1_select. grd -Crainbow -L-5000/3200 -Z > tess. cpt • grdimage -R 10/30/20/50. input_dataEtopo 1_select. grd -Ctess. cpt -JM 5 i -O -K >> Mepo_Es 1_Topo. ps • pscoast -R 10/30/20/50 -JM 5 i -I 1 -Na -W -Dh -P -Ba -O -K >> Mepo_Es 1_Topo. ps • grdcontour. input_dataEtopo 1_select. grd -R 10/30/20/50 -JM 5 i -C 500 -O >> Mepo_Es 1_Topo. ps • c: Programsgsgs 9. 25bingswin 64 -s. DEVICE=jpeg -d. Use. Crop. Box -r 300 s. Output. File=Mepo_Es 1_Topo. jpg Mepo_Es 1_Topo. ps

Indice degli argomenti da trattare durante il laboratorio 1 Installazione e primo utilizzo del pacchetto software GMT (General Mapping Tools). Installare anche database di topografia e coste e fiumi. Creazione di mappe con topografia e coste. 2 Creazione di mappe di: anomalia di gravita’ , Bouguer e topografia per zona a scelta dello studente. Recupero dati dal database ICGEM 3 Modellazione di un corpo tramite un prisma. Analisi delle proprieta’ del campo in funzione delle dimensioni e profondita’ del corpo. Vedi esempi sul testo Hinze, von Freese, Saad, pagina 183. 4 Modellazione di una area tramite tesseroidi 5 Calcolo di anomalia isostatica

Introduzione l’ambiente di sviluppo Mat. Lab • Matlab e’ un ambiente di sviluppo molto potente della Mathworks di larghissimo utilizzo. E’ un sistema integrato per la progettazione, esecuzione e rappresentazione grafica dei risultati di un programma di lavoro. Un programma di lavoro consiste in una serie di comandi predefiniti. Uno script e’ un file contenente i comandi di linea (le routine) che vengono eseguiti in sequenza. • Matlab mette a disposizione una moltitudine di routine classificati a seconda dell’applicativo, come: • - analisi statistica dei dati • -Soluzione di problemi numerici, di inversione, soluzone di equazioni differenziali • -Acquisizione dati in laboratorio • GIS: creare mappa geografica di valori misurati in campagna, scegliere la proiezione cartografica, sovrapposizione di dati di diversa natura: Immagine satellitare, temperatura al suolo a 100 m di profondita’, isolinee della topografia, rete stradale. • Soluzione numerica di problemi ad elementi finiti o alle differenze finite. • Pacchetto Signal Processing. Analisi spettrale, filtraggio di sequenze temporali o nello spazio.

Ambiente integrato Matlab

Cartella di lavoro Create una cartella di lavoro sulla quale raccogliete i vostri script e dati. Per il nome della cartella non utilizzare spazi. Per aggiungere la visibililta’ di un’altra cartella: -> set path

Finestre di Matlab • Current folder: mostra contenuto della cartella aperta. I risultati e le chiamate degli script si riferiscono a questa cartella, se non diversamente specificato • Command window: qui vengono immessi comandi di linea e viene scritto il risultato di ogni elaborazione. Digitando un comando il sistema cerca uno script dello stesso nome nella cartella aperta, oppure in tutte le cartelle incluse nella definizione di path. Il path permette all’interprete di matlab di trovare script di libreria e di sistema. Matlab assegna di default un path valido per lavorare. Le cartelle personali di lavoro devono essere aggiunte manualmente nel path. • L’interprete cerca i files solo nelle cartelle assegnate, vale per gli script e per i dati. Durante l’installazione di matlab vengono definite le cartelle accessibili. Tutti i comandi di Matlab sono salvati in una cartella specifica. Noi definiremo i propri script che salveremo e potremo richiamare nel workspace.

Help • Matlab dispone di un’assistenza ai comandi integrata. Questa funzione e’ fondamentale nell’apprendimento del linguaggio di programmazione. >> help Per informazioni su una funzione specifica: >> help sin Per documentazione piu’ ampia: >> doc sin Per cercare documentazione su una funzione: >> doc e tab sin (appare un menu di ricerca)

Scripts - introduzione • Script: • Una serie di comandi raccolti in un file. I comandi vengono eseguiti consecutivamente • Lo script viene descritto nell’editor di Matlab e salvato come file Matlab, caratterizzato dall’ estensione miofile. m • Da riga di commando: >> edit ciao. m Oppure:

Editor di scripts

Scripts- nota bene • Nello script le righe che iniziano con % vengono interpretate commenti • Le prime righe commentate vengono utilizzate dal “help” per fornire informazioni sul comando rappresentato dallo script • Le variabili create durante l’esecuzione dello script sono visibili nello workspace. Rimangono in memoria anche dopo il termine dell’esecuzione dello script.

Variabili • Iniziamo a lavorare con variabili. Le variabili hanno lo scopo di essere contenitori di numeri con le quali definiamo espressioni matematiche. • Il problema da risolvere trova una formulazione matematica generale senza la necessita‘ di assegnare esplicitamento un numero alla variabile. • Al momento dell‘esecuzione dello script, dobbiamo assegnare un valore esplicito ad una variabile prima di utilizzarla. Altre variabili vengono assegnate al momento del calcolo. • Tipi di variabili: numeri (64 -bit double), stringhe di caratteri (16 -bit char), inoltre variabili simboliche, complesse, integer. • Comunemente utilizzeremo array di numeri o caratteri

Nome delle Variabili • • • Per assegnare un valore ad una variabile: >> var 1= 2. 7183 >> stringa=‘franz’ Nomi delle variabili: il primo carattere deve essere una lettera, poi combinazione di lettere e numeri. Case-sensitive: distinzione di minuscole e maiuscole! Variabili di sistema, da non utilizzare: i e j potrebbere indicare numeri complessi pi greco: pi= 3. 1416 ans: ultimo valore Inf e –Inf: positivo e negativo infinito Na. N: ‘not a number’ nei calcoli viene trattato come numero inesistente.

Array • Molto spesso necessitiamo di raccogliere un insieme di numeri in una unica variabile. Definiamo un array di numeri. • Esempio: tutti i valori della temperatura misurata nell’arco di un giorno. • In Matlab possiamo anche definire un array di stringhe e/o numeri, ed e’ un cell array. • Esempio: tutti i nomi degli studenti del corso

Vettore riga • >> riga = [1 3 5 6. 3 7. 5 3. 2] • • >> riga = [1 3 5 6. 3 7. 5 3. 2] Output in Command window: riga = 1. 0000 3. 0000 5. 0000 6. 3000 7. 5000 3. 2000 • Workspace:

Vettore colonna • >> col = [2; 4; 6; 8] • • • >> Output in Command window: col = 2 4 6 8

Operazioni sugli elementi di un array • Le funzioni che operano su scalari, solitamente operano anche su array. • >> f=[0 pi/4 pi/2 3*pi/4 ] • >> r=sin(f) • Equivale a: • >> r = [sin(0) sin(pi/4) sin(pi/2) sin(3*pi/4) ] • Invece le operazioni (* / ^) distinguono fra operazioni su singolo elemento dell’array o operazioni fra array.

Operatori aritmetici su singolo elemento • Le operazioni (* / ^) applicate su un array effttuano operazioni di calcolo matriciale. Se invece intendiamo operazione elemento per elemento, e’ necessario anteporre un punto:

Grafico di una sequenza • Un grafico di una serie di valori, come quelli della temperatura richiede una serie di coppie di valori, x e y. >>x=linspace(0, 20); >>y=3*sin(2*pi/ 10 *x); >>plot(x, y); >>title('Sinusoide'); >>xlabel('x(m)’; >>ylabel('y(m)');

L’anomalia di gravità di una sfera- Matlab %% calculate the gravity effect of a sphere with radius R at depth D and located at x=0 along a profile x=-100: 0. 5: 100; % abscissa of the profile G=6. 67*10^-11; % Gravit. Const ms 2_m. Gal=10^5; % conversion m/s^2 2 m. Gal dens=1000; % density of the body kg/m 3 R=20; % radius of the sphere [m] D=20; % depth of the center of the sphere [m] x 0=0; % coordinate in the x direction of the center of sphere [m] gz=(4/3*pi*G*ms 2_m. Gal*dens*R^3). /(3*D^2. *(1+((x-x 0). ^2/D^2)). ^(3/2)); figure plot(x, gz) % command for plotting a function of 1 variable title('gravity anomaly due to a sphere') xlabel('x [m]') ylabel('gravity anomaly [m. Gal]')

Esercizi • Variare i parametri della sfera (raggio, profondità, densità). Comparare i profili delle anomalie. Cosa si osserva? • Calcolare anomalie dovute ad altri corpi vedi tabella 3. 2 (pag. 49 Hinze et al. ). Cilindro orizzontale infinito e inclined sheet (pag 50 tab. 3. 3).

Indice degli argomenti da trattare durante il laboratorio 1 Installazione e primo utilizzo del pacchetto software GMT (General Mapping Tools). Installare anche database di topografia e coste e fiumi. Creazione di mappe con topografia e coste. 2 Creazione di mappe di: anomalia di gravita’ , Bouguer e topografia per zona a scelta dello studente. Recupero dati dal database ICGEM 3 Modellazione di un corpo tramite un prisma. Analisi delle proprieta’ del campo in funzione delle dimensioni e profondita’ del corpo. Vedi esempi sul testo Hinze, von Freese, Saad, pagina 183. 4 Calcolo Radice Isostatica 5 Modellazione della radice isostatica tramite tesseroidi e calcolo dell’ anomalia isostatica

Contenuti della lezione • Topografia, Anomalia di gravità e Anomalia di Bouguer • Calcolo della radice isostatica con GMT • Plot della radice con GMT

Tesseroids: il programma • Scaricate dal sito il file tesseroids-1. 2. 1. zip e scompattatelo in C: MEPO • https: //github. com/leouieda/tesseroids/releases/tag/v 1. 2. 1 • Tesseroids è una raccolta di programmi scritti da Leonardo Uieda (Observatorio Nacional di Rio) • Consentono di calcolare Potenziale e derivate prime e seconde dovute ad un modello di tesseroidi

Create le carte di Topografia, Anomalia di Gravità e Bouguer per un’area orogenica: • Topografia sviluppata fino a grado e ordine 280 • Gravity_anomaly_sa fino a g/o 280 calcolata a 4000 m di quota • Gravity_anomaly_bg fino a g/o 280 calcolata a 4000 m di quota • Posizionate i files scaricati nella cartella dati_input • Lo script sul quale lavorate xxx. bat deve essere posizionato nella cartella tesseroids-1. 2. 1. La cartella con i dati “dati_input” deve essere una sotto-cartella della cartella tesseroids-1. 2. 1.

Convertire i file da tabella in griglia. Nell’esempio qui sotto il header dei singoli file e’ stato eliminato. • : : convert xyz datafile to gmt grid. Calc • xyz 2 grd dati_inputTopo 280 NOHEAD. xyz –Gdati_inputtopo. grd -I 0. 2 -R 10/15/45/50 -V • xyz 2 grd dati_inputFAgoco 05 s-NOHEAD. xyz -Gdati_inputFA. grd -I 0. 2 -R 10/15/45/50 -V • xyz 2 grd dati_inputBGgoco 05 s-NOHEAD. xyz -Gdati_inputBG. grd -I 0. 2 -R 10/15/45/50 -V

Creare le tre mappe dei campi • • : : plot the three maps psbasemap -R 10/15/45/50 -JM 5 i -P -Ba -B+t"Topo" -K > Topo. ps grd 2 cpt dati_inputtopo. grd -Z -L-300/2500 -Crainbow > tess. cpt grdimage -R 10/15/45/50 -JM 5 i dati_inputtopo. grd -Ctess. cpt -O -K >> Topo. ps • psscale -Dx 6 i/1 i+w 4 i/0. 5 i -Ba -Ctess. cpt -O -K >> Topo. ps • pscoast -R 10/15/45/50 -JM 5 i -I 1 -P -Na -W -O >> Topo. ps • • set grid=dati_inputFA. grd set out=FA. ps psbasemap -R 10/15/45/50 -JM 5 i -P -Ba -B+t"FA" -K > %out% grd 2 cpt %grid% -Z -L-100/100 -Crainbow > tess. cpt grdimage -R 10/15/45/50 -JM 5 i %grid% -Ctess. cpt -O -K >> %out% psscale -Dx 6 i/1 i+w 4 i/0. 5 i -Ba -Ctess. cpt -O -K >> %out% pscoast -R 10/15/45/50 -JM 5 i -I 1 -P -Na -W -O >> %out%. . set grid=dati_inputBG. grd set out=BG. ps psbasemap -R 10/15/45/50 -JM 5 i -P -Ba -B+t"BG" -K > %out% grd 2 cpt %grid% -Z -L-200/20 -Crainbow > tess. cpt grdimage -R 10/15/45/50 -JM 5 i %grid% -Ctess. cpt -O -K >> %out% psscale -Dx 6 i/1 i+w 4 i/0. 5 i -Ba -Ctess. cpt -O -K >> %out% pscoast -R 10/15/45/50 -JM 5 i -I 1 -P -Na -W -O >> %out%

Topografia e anomalie di gravità Topografia Anomalia di Gravità Bacino molasse Alpi Pianura e adriatico

Topografia- Bouguer • Topografia Bouguer Bacino molasse Alpi Pianura e adriatico

Carta delle anomalie di Bouguer • Dominata in termini di ampiezza da anomalie sotto le Alpi • Sotto le Alpi c’è una mancanza di massa • Vogliamo verificare se è ‘giustificabile’ dall’effetto isostatico

Isostasia secondo Airy Bilancio isostatico a terra ed esplicitazione della radice b 1 in funzione degli altri parametri Da Wikipedia Bilancio isostatico a mare ed esplicitazione della radice b 2 in funzione degli altri parametri

Il calcolo della radice isostatica in GMT • Set input=dati_inputtopo. grd • Set flex=dati_inputflex. grd • Sopra abbiamo definito due variabili che possono essere utilizzate nei comandi gmt. grdflexure %input% -D 3200/2670/1035 -E 0 k -G%flex% -N+l -Z 30000 –fg %input% nome del file di topografia -D 3300/2700/1035 densita’ mantello, crosta, acqua -E 0 k spessore elastico in km. Per isostasia Airy: Te=0 km -G%flex% nome del grid di flessione. -N+l non modificare il file in input (non togliere media, piano inclinato o altro) -Z 30000 profondita’ in metri del livello Moho di flessione senza carico (Profondita’ di riferimento) • –fg lavoriamo in coordinate geografiche (e non metriche). • •

Il calcolo della radice isostatica secondo Airy in GMT (valido solo per zone emerse) • Set input=dati_inputtopo. grd • Set flex=dati_inputflex. grd • grdmath %input% -2670 MUL 530 DIV -30000 ADD = %flex% • %input% nome del file di topografia • -2670 densita’ crosta, 530 contrasto densita’ mantello/crosta • %flex% nome del grid di flessione. • -30000 profondita’ in metri del livello Moho di flessione senza carico (Profondita’ di riferimento) • ADD somma. MUL moltiplica. DIV divide-

La radice di Airy nelle Alpi Orientali

Indice degli argomenti da trattare durante il laboratorio 1 Installazione e primo utilizzo del pacchetto software GMT (General Mapping Tools). Installare anche database di topografia e coste e fiumi. Creazione di mappe con topografia e coste. 2 Creazione di mappe di: anomalia di gravita’ , Bouguer e topografia per zona a scelta dello studente. Recupero dati dal database ICGEM 3 Modellazione di un corpo tramite un prisma. Analisi delle proprieta’ del campo in funzione delle dimensioni e profondita’ del corpo. Vedi esempi sul testo Hinze, von Freese, Saad, pagina 183. 4 Calcolo Radice Isostatica 5 Modellazione della radice isostatica tramite tesseroidi e calcolo dell’ anomalia isostatica

Contenuti del modulo • Calcolo dell’effetto gravimetrico della radice isostatica (modellazione diretta) • Confronto con il dato osservato della Carta di Bouguer • La carta dei residui

Modellazione diretta - 1 • L’effetto di gravità prodotto da alcune geometrie elementari è noto. • Ad esempio conosciamo l’effetto di gravità di un cilindro, di una sfera, di un prisma • Per il campo di gravità come altri campi (elettrostatico, magnetico) vale il principio di sovrapposizione degli effetti: il campo dovuto ad una serie di masse è uguale alla somma dei singoli campi dovuti a ciascuna massa • Da questa considerazione si configura la possibilità di calcolare il campo dovuto a complesse geometrie attraverso la somma degli effetti di geometrie elementari.

Modellazione diretta - 2 • Modellazione diretta: processo con il quale si cerca di riprodurre il campo osservato mediante il campo calcolato a partire da un modello di densità da noi ipotizzato. Zhang et al. , 2014

Effetto gravimetrico della radice – modello di riferimento • Possiamo calcolare l’effetto sul campo di gravità delle variazioni laterali di massa dovute a questa radice Modello di Riferimento Crosta 2670 Modello di cui vogliamo calcolare l’effetto su g Crosta 2670 Mantello 3200 Modello discretizzato a tesseroidi La differenza tra gli effetti gravimetrici di questi due modelli definisce un’anomalia di gravità che possiamo confrontare con le anomalie osservate nella carta di Bouguer

Effetto gravimetrico della radicerappresentazione contrasto di densità - gz = gz 0 m. Gal x Crosta 2670 Mantello 3200 = -530

Tesseroidi: prismi in coordinate sferiche • «Building block» per costruire modelli complessi • Un tesseroide è delimitato da 2 piani in longitudine , 2 in latitudine e 2 superfici sferiche • A ciascun tesseroide è assegnata una densità

Tesseroids: il programma N righe N tesseroidi • Tesseroids è una raccolta di programmi scritti da Leonardo Uieda (Observatorio Nacional di Rio) • Consentono di calcolare Potenziale e derivate prime e seconde dovute ad un modello di tesseroidi • Un modello è definito da un file testuale così strutturato: Lon 1 Lon 2 Lat 1 Lat 2 top bottom contrasto_densità 9. 9 10. 1 49. 9 50. 1 -30000 -31425 -530 10. 1 10. 3 49. 9 50. 1 -30000 -31032 -530. . .

Tesseroids: il programma • Il programma non ha interfaccia grafica, ma viene lanciato da linea di comando • Come gia’ visto in precedenza, e’ necessario creare un file. bat che contiene lo script per l’esecuzione dei comandi in sequenza.

Tesseroids – Script di calcolo - Create un file calcolo. bat (in unix calcolo. sh) - Vengono chiamati tre programmi in sequenza e sono tessgrd e tessgz - Ciascun programma ha degli input passati sia come file che come parametri (es: –r-5/5/-5/5. . . ) - Gli output di ciascun programma sono dopo il >. - Salvatelo nella cartella tesseroids.

Tesseroids – Script di calcolo- seconda riga •

Tesseroids – Script di calcolo – terza riga bintessgz. exe dati_inputmodello_densita. txt < dati_inputgrid. Calc. Test. txt > dati_inputeffetto. Test. txt bintessgz. exe chiama tessgz programma che calcola l’effetto gravimetrico del modello. In input accetta due file il modello tesseroidi elabmodello_densita. txt e il grid di calcolo dati_inputgrid. Calc. Test. txt Contenuto del file modello_densita. txt da creare nella cartella dati_input -1 1 0 -5000 -300

Cartella di lavoro tesseroids 1. 2. 1 Cartella di lavoro ‘dati_input’ con il file modello_densita. txt Cartella di Tesseroids con la cartella di lavoro ‘dati_input’ e lo script calcolo. bat

Il file grid. Calc. Test. txt

Output da tesseroids- effetto gravimetrico # gz component calculated with tessgz 1. 2. 0: # local time: Wed Dec 06 11: 20: 53 2017 # model file: modello_densita. txt (1 tesseroids) # GLQ order: 2 lon / 2 lat / 2 r # Use recursive division of tesseroids: True # Distance-size ratio for recusive division: 1. 5 # Grid generated with tessgrd 1. 2. 0: # local time: Wed Dec 06 11: 20: 53 2017 # args: -r-5/5/-5/5 -b 26/26 -z 4000 # grid spacing: 0. 400000 lon / 0. 400000 lat # total 676 points -5 -5 4000 -0. 0564012886092919 -4. 6 -5 4000 -0. 0593200806498219 -4. 2 -5 4000 -0. 0624165460154603 -3. 8 -5 4000 -0. 065673518488156 -3. 4 -5 4000 -0. 069058821950435 -3 -5 4000 -0. 0725207791573518 -2. 6 -5 4000 -0. 0759841250443243 -2. 2 -5 4000 -0. 0793476511357376 -1. 8 -5 4000 -0. 0824853632481009 -1. 4 -5 4000 -0. 085252914150674 -1 -5 4000 -0. 0875001671208956 -0. 6 -5 4000 -0. 089088889144073 Header File a 4 colonne con longitudine, latitudine, quota di calcolo e valore del funzionale calcolato

Converti tabella da tesseroids in griglia • gmtconvert dati_inputeffetto. Test. txt i 0, 1, 3 > tmp Seleziona prime due e la terza colonna • xyz 2 grd tmp -Gdati_inputeffetto. grd I 0. 4/0. 4 -R-5/5/-5/5 -V • grdinfo dati_inputeffetto. grd # gz component calculated with tessgz 1. 2. 0: # local time: Mon Jan 08 21: 50 2018 # model file: elabmodello_densita. txt (1 tesseroids) # GLQ order: 2 lon / 2 lat / 2 r # Use recursive division of tesseroids: True # Distance-size ratio for recusive division: 1. 5 # Grid generated with tessgrd 1. 2. 0: # local time: Mon Jan 08 21: 50: 49 2018 # args: -r-5/5/-5/5 -b 26/26 -z 4000 # grid spacing: 0. 400000 lon / 0. 400000 lat # total 676 points -5 -5 4000 -0. 0564012886092919 -4. 6 -5 4000 -0. 0593200806498219 -4. 2 -5 4000 -0. 0624165460154603

Sovrapporre basemap con coste all’effetto tesseroidi (mappa a colori) • psbasemap -R-5/10/-5/5 -JM 6 i -P -Ba -B+t"Test Tesseroide" -K > Test_tess. ps • grd 2 cpt dati_inputeffetto. grd -Crainbow -Z > tess. cpt • grdimage dati_inputeffetto. grd -R-5/10/-5/5 -Ctess. cpt -JM 6 i -O -K >> Test_tess. ps • grdcontour dati_inputeffetto. grd -R-5/10/-5/5 -JM 6 i -C 10 -O -K >> Test_tess. ps • psscale -Dx 6. 5 i/0 i+w 4 i/0. 1 i -Ba -Ctess. cpt -O -K >> Test_tess. ps • pscoast -R-5/10/-5/5 -JM 6 i -I 1 -P -Ba -O -W >> Test_tess. ps

Effetto gravimetrico di un tesseroide

Tesseroids applicato alle Alpi • Come passare da un grid XYZ a un modello a tesseroidi? • Fra i programmi disponibili in Tesseroids 1. 2. 1 c’è un convertitore da file grid ascii xyz a tesseroidi. • In pratica attorno ad ogni punto del grid viene costruito un tesseroide di grandezza uguale al passo • Dobbiamo solo impostare il contrasto di densità e la quota di riferimento ed i tesseroidi verranno di fatto definiti attorno alla superficie di riferimento in questo modo: Oscillazione Moho +530 -530 Profondità di Riferimento 30000 m

Tesseroids – Script di calcolo : : tesseroid calculation set flex_xyz=dati_inputflex. txt grd 2 xyz %flex% > %flex_xyz% set tess=dati_inputtessflex. txt set points=dati_inputgrid. Calc. txt set gflex=dati_inputflex. Effect. txt bintessmodgen. exe -s 0. 1/0. 1 -d 530 -z-30000 -v < %flex_xyz% > %tess% bintessgrd. exe -r 10/15/45/50 -b 26/26 -z 4000 -v > %points% bintessgz. exe %tess% < %points% > %gflex% - Vengono chiamati tre programmi in sequenza e sono tessmodgen, tessgrd e tessgz - Ciascun programma ha degli input passati sia come file che come parametri (es: -s 0. 1/0. 1 –d 530. . . ) - Gli output di ciascun programma sono dopo il >.

Tesseroids – Script di calcolo- prima riga set flex=dati_inputflex. grd set flex_xyz=dati_inputflex. txt set tess=dati_inputtessflex. txt grd 2 xyz %flex% > %flex_xyz% bintessmodgen. exe -s 0. 1/0. 1 -d 530 -z-30000 -v < %flex_xyz% > %tess% bintessmodgen. exe = chiamata del programma che converte file. xyz in modello tesseroidi. -s 0. 1/0. 1= parametro in cui sono indiacate la spaziature del grid. xyz nelle direzioni longitudine e latitudine -d 530= contrasto di densità della superficie di cui stiamo calcolando l’effetto gravimetrico (crosta/mantello) -z-30000 profondità di riferimento Attraverso la sintassi < %flex_xyz% diciamo che l’input è il file flex. txt Attraverso la sintassi > %tess% diciamo che salvi le elaborazioni nel file dati_inputtessflex. txt

Tesseroids – Script di calcolo- seconda riga set points=dati_inputgrid. Calc. txt bintessgrd. exe -r 10/15/45/50 -b 51/51 -z 4000 -v > %points% bintessgrd. exe chiamiamo il programma tess grd che consente di definire la griglia di calcolo: ovvero dove e a che quota vogliamo calcolare gli effetti del nostro modello. -r 10/15/45/50 definisce gli estremi della grid di calcolo: da 10 a 15 in longitudine e da 45 a 50 in latitudine -b 51/51 numero di divisioni del grid lungo la longitudine e latitudine; si può ricavare la spaziatura da: (15 -10)/(b-1)=0. 1° lungo la longitudine -z 4000 è la quota di calcolo -v indica la modalità verbose del programma, così gli eventuali errori vengono stampati nella linea di comando L’output viene scritto nella directory dati_input, con nome file di grid. Calc. txt

Tesseroids – Script di calcolo – terza riga set tess=dati_inputtessflex. txt set points=dati_inputgrid. Calc. txt set gflex=dati_inputflex. Effect. txt bintessgz. exe %tess% < %points% > %gflex% bintessgz. exe chiama tessgz programma che calcola l’effetto gravimetrico del modello. In input accetta due file il modello tesseroidi dati_inputtessflex. txt e il grid di calcolo dati_inputgrid. Calc. txt

Output da tesseroids- effetto gravimetrico della Moho di Airy # gz component calculated with tessgz 1. 2. 0: # local time: Sun Dec 04 15: 35: 55 2016 # model file: tess. Airy. txt (676 tesseroids) # GLQ order: 2 lon / 2 lat / 2 r # Use recursive division of tesseroids: True # Distance-size ratio for recusive division: 1. 5 # Grid generated with tessgrd 1. 2. 0: # local time: Sun Dec 04 15: 35: 55 2016 # args: -r 10/15/45/50 -b 26/26 -z 4000 # grid spacing: 0. 200000 lon / 0. 200000 lat # total 676 points 10 45 4000 -10. 9960136985812 10. 2 45 4000 -11. 2388002259674 10. 4 45 4000 -10. 7237737154379 10. 6 45 4000 -10. 1463749762208 10. 8 45 4000 -9. 92683017256846 11 45 4000 -10. 1089202371021 11. 2 45 4000 -10. 4852420813538 11. 4 45 4000 -10. 7631533654661 11. 6 45 4000 -10. 7126752078471 11. 8 45 4000 -10. 2591541451825 12 45 4000 -9. 50551013353708 12. 2 45 4000 -8. 68914331793644 Header File a 4 colonne con longitudine, latitudine, quota di calcolo e valore del funzionale calcolato. Assomiglia al file. xyz

Trasformare la tabella in formato griglia leggibile da GMT : : convert xyz datafile to gmt grid. Calc gmtconvert %gflex% -i 0, 1, 3 > tmp : : Seleziona prime due e la terza colonna set gflexgrid=dati_inputflex. Effect. grd xyz 2 grd tmp -G%gflexgrid% -I 0. 1 -R 10/15/45/50 -V Grid step del grid Estensione areale %gflexgrid% del grid %gflexgrid%

Effetto della radice isostatica

Carta dei Resiudi • Il residuo sarà semplicemente data dalla differenza tra il grid Bouguer osservato e il grid effetto Airy. • ATTENZIONE: In questo caso abbiamo calcolato l’anomalia della radice su un area più piccola rispetto alla carta delle anomalie di Bouguer scaricata da ICGEM. Quindi tagliamo quest’ultima con gli estremi dell’area desiderata con il comando grdcut • set iso=dati_inputIso_res. grd • grdcut dati_inputGrav. Anom_BG. grd -Gdati_inputBG_grav_cut. grd -R 10/15/45/50 -V • grdmath -R 10/15/45/50 dati_inputBG_grav_cut. grd %gflexgrid% SUB = %iso% • : : graph of gravity residual • set out=Isores. ps • psbasemap -R 10/15/45/50 -JM 5 i -P -Ba -B+t"Iso Residual" -K > %out% • grd 2 cpt %iso% -Z -L-100/100 -Crainbow > tess. cpt • grdimage -R 10/15/45/50 -JM 5 i %iso% -Ctess. cpt -O -K >> %out% • psscale -Dx 6 i/1 i+w 4 i/0. 5 i -Ba -Ctess. cpt -O -K >> %out% • pscoast -R 10/15/45/50 -JM 5 i -I 1 -P -Na -W -O >> %out%

Residuo Isostatico