Geotermometro clino/orto-pirosseno¶
Vedremo qui di mettere in pratica la teoria discussa nelle lezioni di termodinamica (rif. lezione termodinamica_17), con l'esempio di un ipotetico geotermometro basato sulla ripartizione di ferro e magnesio tra due pirosseni, di cui uno monoclino e l'altro ortorombico.
Per chi volesse replicare i dati usati per la calibrazione del geotermometro, questi sono stati simulati usando il programma Perplex, con il database termodinamico hp622ver.dat per gli end-member, i modelli di soluzione Opx(JH) e Cpx(JH), e per una composizione globale (in moli) pari a:
CaO 0.4; MgO 1.4; FeO 0.2 e SiO2 2.0.
I dati ottenuti dalla simulazione con Perplex (frazioni molari di Fe nel clinopirosseno, nell'ortopirosseno e la temperatura corrispondente) sono contenuti nel file Excel geotermometro.xlsx che deve essere presente nel folder da cui viene lanciato il programma o questo notebook.
Il notebook funziona sulla base del programma geotermometro.py che lanceremo nelle celle successive.
import inspect as ins
%matplotlib inline
%run geotermometro.py
Dopo aver lanciato il programma (cella precedente) devono essere caricati i dati per la calibrazione del geotermometro, usando la funzione read_data, il cui codice è:
print(ins.getsource(read_data))
Dal punto di vista dell'attività di programmazione, l'interesse per questa funzione sta nell'uso delle funzioni della libreria pandas per leggere dataset contenuti in file di vario formato (tra cui Excel) e di manipolarli nel modo opportuno. In questo caso si usa la funzione read_excel della libreria pandas per caricare il file ed assegnare le variabili t_list (T: temperatura), c_list (CPX: frazione molare del ferro nel clinopirosseno) e o_list (OPX: frazione molare del ferro nell'ortopirosseno) con i dati contenuti nel file.
I dati caricati devono essere condivisi tra funzioni diverse del programma e, per questo motivo, vengono salvati entro una variabile dt che è l'istanza di una semplice classe data il cui codice è
class data():
def __init__(self):
self.t_list=np.array([])
self.c_list=np.array([])
self.o_list=np.array([])
self.reg=np.array([])
self.data_flag=False
self.reg_flag=False
self.mnr=0.
self.mxr=0.
self.ll_dat=0
def set_reg(self,reg,mnr,mxr):
self.reg_flag=True
self.reg=reg
self.mnr=mnr
self.mxr=mxrLa variabile dt è poi creata con l'assegnazione
dt=data()
In questo modo, la lista dei valori di temperatura viene ad esempio conservata nella variabile dt.t_list. La classe prevede anche un metodo (set_reg) per assegnare i parametri della regressione lineare del rapporto CPX/OPX in funzione della temperatura (che sono usati in varie funzioni del programma)
read_data()
La lista dei dati di calibrazione viene stampata come illustrato nella cella precedente da una funzione print_data chiamata da read_data, il cui codice è:
print(ins.getsource(print_data))
La funzione print_data sfrutta alcune funzionalità della libreria pandas per strutturare velocemente delle tabelle di dati per una chiara visualizzazione. Per esempio, supponiamo di creare una lista di valori di una variabile indipendente x e del suo quadrato (y=x^2); possiamo costruire una tabella con il codice che segue:
n_val=5
x_ini=0.
x_fin=10.
# Liste x e y
x_list=np.linspace(x_ini,x_fin,n_val)
y_list=x_list**2
tab=(x_list,y_list)
data=pd.DataFrame(tab, index=['X', 'Y'])
print(data)
Volendo scrivere la tabella in due colonne, trasponiamo il DataFrame data scrivendo semplicemente:
data=data.T
print(data)
Per eliminare i numeri di linea un modo semplice è:
print(data.to_string(index=False))
che sfrutta il metodo to_string per convertire i dati numerici della tabella in righe di tipo carattere, rimuovendo l'indice (index=False che corrisponde al numero di linea). Esistono poi molti attributi e metodi scritti per i DataFrame di pandas che consentono di operare sui dati. Per esempio i due seguenti:
data['Y'].values # accede ai dati della colonna "Y"
data['X'].mean() # metodo per calcolare la media dei dati di colonna 'X'
Tornando al geotermometro, la calibrazione viene fatta dalla funzione calibrazione: questa calcola il rapporto (ratio) tra le frazioni molari del ferro nel clinopirosseno e nell'ortopirosseno per ogni valore di temperatura (dt.t_list) e poi fa un fit lineare di T in funzione di ratio: T=m\*ratio + q. La funzione polyfit della libreria Numpy esegue il fit; i risultati (m e q) vengono passati alla variabile reg che viene poi salvata nella variabile dt (istanza della classe data), usando il metodo set_reg (in modo che le altre funzioni possano accedervi; l'uso di questa tecnica evita il ricorso alle dichiarazioni global di variabili create all'interno di funzioni e poi utilizzate all'esterno di esse; come già accennato in altri notebook, l'uso di dichiarative global è sconsigliato). date un'occhiata alla funzione calibrazione:
print(ins.getsource(calibrazione))
Oltre che fare il fit, la funzione si preoccupa di mostrare un grafico con i dati originali e la retta di regressione. Lanciamola:
calibrazione()
Come detto nella lezione termodinamica_17, con riferimento alla reazione di scambio:
$${\rm Mg}^c + {\rm Fe}^o \leftrightarrow {\rm Mg}^o + {\rm Fe}^c$$(il simbolo $A^\alpha$ indica la specie chimica $A$ nella fase $\alpha$; $c$ sta per clinopirosseno e $o$ sta per ortopirosseno) possiamo scrivere la costante di equilibrio (costante a temperatura e pressione costante) come
$$K = \frac{x_{\rm Mg}^{o} x_{\rm Fe}^{c}}{x_{\rm Mg}^{c} x_{\rm Fe}^{o}}= \frac{x_{\rm Fe}^{c}/x_{\rm Fe}^{o}}{x_{\rm Mg}^{c}/x_{\rm Mg}^{o}}$$Come visto nelle lezioni, esiste una relazione tra $K$ e la temperatura:
$$\log K= -\frac{\Delta\mu^0}{RT}$$dove $\Delta\mu^0$ è la variazione di energia libera molare associata alla reazione medesima e dipende dalle energie libere molari dei termini puri alla temperatura e pressione determinata. Essendo una variazione di energia libera a P e T costanti, sappiamo che questa può essere scritta come:
$$\Delta\mu^0 = \Delta h^0 - T\Delta s^0$$e quindi:
Facendo un fit lineare di $\log K$ in funzione di $1/T$ otteniamo allora l'entalpia della reazione di scambio (coefficiente angolare $m=-\Delta h^0/R$) e l'entropia (intercetta $q=\Delta s^0/R$). La dipendenza di $\log K$ da $1/T$ è lineare ammesso che l'entalpia e l'entropia di reazione non varino apprezzabilmente con la temperatura.
Per chi sia interessato, ecco il codice che calcola l'entalpia e l'entropia della reazione di scambio a partire dai dati usati per la calibrazione del geotermometro (la funzione è relativamente complicata da istruzioni che eseguono un sampling conveniente dei dati di temperatura, essenzialmente per motivi grafici; sorvoliamo su questi aspetti):
print(ins.getsource(log_ratio))
Ed eccone il semplice uso:
log_ratio()
Siamo adesso pronti a utilizzare il geotermometro su ipotetici dati: abbiamo per esempio misurato la frazione molare del ferro in un clinopirosseno e in un ortopirosseno che sono a contatto in un dato campione di roccia, e vogliamo avere una stima della temperatura a cui questa si è equilibrata. Usiamo allo scopo la funzione temperatura, passandogli come argomenti i valori delle frazioni molari misurate.
Il codice della funzione è molto semplice: usa i parametri della retta di regressione ottimizzati dalla funzione calibrazione e calcola direttamente la temperatura note le frazioni molari. Esegue anche un plot per vedere come si colloca il nostro dato $\left(x_{\rm Fe}^{c}, x_{\rm Fe}^{o}\right)$ rispetto ai dati usati per la calibrazione:
print(ins.getsource(temperatura))
Usiamo la funzione:
temperatura(0.050,0.134)
Un'idea dell'errore associato a questa stima la abbiamo guardando all'errore medio del fit dei dati di calibrazione, calcolato con la funzione error:
print(ins.getsource(error))
error()
Possiamo quindi aspettarci un errore di circa 4 K, dovuto al fit dei dati di calibrazione. Poi può esserci un errore legato all'incertezza delle nostre misure sul campione specifico. Per questo è improbabile (e sconsigliabile) avere una sola misura delle concentrazioni della nostra specie chimica nelle due fasi a contatto... E' normale fare misurazioni ripetute in punti adiacenti del campione per poi mediare su tutte le stime di temperatura ottenuta. Supponiamo di aver fatto più misure e di averle riportate in due colonne di un file Excel user_data.xlsx (presente nel folder di questo notebook); usiamo la funzione temp_serie per stimare le corrispondenti temperature e mediare i risultati.
Il codice della funzione è:
print(ins.getsource(temp_serie))
Anche questa funzione usa le funzioni della libreria pandas per leggere i dati del nostro file e per costruire una tabella che riassuma dati e risultati:
temp_serie()
Con quei dati abbiamo una stima della temperatura intorno agli 830K, con un'incertezza di 9K dovuta alla dispersione dei dati medesimi. Questa incertezza va combinata con quella relativa al fit di calibrazione (4K)...
np.sqrt(9**2+4**2).round(1)
Quindi siamo intorno ai 10K di possibile errore sulla stima della temperatura. Ma c'è dell'altro... I nostri dati (asterischi in rosso nel grafico sopra) indicano che forse il nostro sistema non ha raggiunto il pieno equilibrio, oppure che la calibrazione è stata fatta su un campione di roccia chimicamente troppo diverso da quello sul quale stiamo lavorando noi... In questo genere di lavori, mai fidarsi troppo di dati che provengano da un'unica tecnica: li si può accettare provvisoriamente in mancanza di altre evidenze e supporti incrociati, ma bisogna sempre essere pronti a una loro revisione.
Lavoriamo di fantasia...¶
Le temperature per le quali abbiamo fatto la calibrazione del geormometro sono conservate nella lista dt.t_list:
dt.t_list
Selezioniamo le ultime 10, e inseriamole nella lista new_t (la sintassi è un po' strana... ma sorvoliamo qui sugli aspetti informatici), supponendo di aver fatto misure per la calibrazione del geotermometro a partire da circa 1060K fino a 1200K e non da 500K come era per i dati originali.
new_t=dt.t_list[-10::1]
print(new_t)
Facciamo lo stesso con i dati di frazione molare, con i quali ci calcoliamo la costante K
new_c=dt.c_list[-10::1] # frazioni molari del ferro nel clinopirosseno
new_o=dt.o_list[-10::1] # frazioni molari del ferro nell'ortopirosseno
new_mg_c=1-new_c # Mg nel clinopirosseno
new_mg_o=1-new_o # Mg nell'ortopirosseno
K=(new_c/new_o)/(new_mg_c/new_mg_o) # costante di equilibrio
print(K.round(3))
Facciamo un plot di della costante di equilibrio in funzione della temperatura:
plt.figure()
plt.plot(K,new_t,"k*")
plt.ylabel("T (K)")
plt.xlabel("K")
plt.show()
I punti ottenuti appaiono abbastanza ben allineati su una retta... Quindi, non conoscendo nulla di termodinamica, perchè non farci una bella regressione lineare?
Pronti!
lin=np.polyfit(K,new_t,1) # regressione lineare sui dati T(K)
k2=np.array([0.40,0.52]) # valori di K e T dal fit, da usarsi nella grafica
t2=np.polyval(lin,k2)
# grafica per il controllo del fit
plt.figure()
plt.plot(K,new_t,"k*",label="dati attuali")
plt.plot(k2,t2,label="fit lineare")
plt.ylabel("T (K)")
plt.xlabel("K")
plt.legend(frameon=False)
plt.show()
Bene... adesso siamo pronti per l'estrapolazione alle temperature più basse, usando il fit lineare che ci siamo costruiti. Dal laboratorio accanto, il nostro collega ci passa i risultati delle sue misure su un campione incognito: le frazioni molari del ferro nel clinopirosseno e nell'ortopirosseno sono rispettivamente 0.010 e 0.162, e ci chiede una stima della temperatura di equilibrio:
cfe=0.010 # dati iniziali
ofe=0.162
cmg=1-cfe
omg=1-ofe
rmg=cmg/omg
rfe=cfe/ofe
k_sample=rfe/rmg
T_est=np.polyval(lin,k_sample) # stima della temperatura sulla base del K misurato
# (k_sample) e dei parametri della retta di regressione
# ottimizzati precedentemente (lin)
print("Stima della temperatura: %5.2f (K)\n" % T_est)
k2=np.array([0.0, 0.55]) # calcolo T per due i estremi (min e max di K)
t2=np.polyval(lin,k2) # sulla base del fit lineare
plt.figure()
plt.plot(K,new_t,"k*",markersize=4,label="dati calibrazione")
plt.plot(k2,t2,label="fit lineare")
plt.scatter(k_sample,T_est,c='red', label="dato incognito")
plt.ylabel("T (K)")
plt.xlabel("K")
plt.legend(frameon=False)
plt.show()
Con il nostro fit lineare, del tipo:
$$T=a+bK$$dove $a$ e $b$ sono le due variabili di regressione ottimizzate sulla base del nostro (limitato) set di dati, abbiamo stimato una temperatura di equilibrio di circa 532K. Adesso però, qualcuno che conosce la termodinamica meglio di noi, ci dice che quel fit lineare che secondo noi, sulla base del grafico, legherebbe $T$ e $K$, non ha alcun fondamento nonostante le apparenze. Quello che dovremmo invece fare, volendo ancora usare un fit lineare, è legare il logaritmo di $K$ al reciproco della temperatura:
$$\frac{1}{T} = a + b\log K$$con $a$ e $b$ da ottimizzarsi sulla base dei nostri dati... Un po' scettici circa questo avvertimento, andiamo a fare questo fit e poi utilizziamolo per rifare la nostra stima della temperatura:
log_new_k=np.log(K) # logaritmo di K
rt_new=1./new_t # reciproco della temperatura
termo_fit=np.polyfit(log_new_k,rt_new,1) # fit di 1/T in funzione di log(K)
log_k_sample=np.log(k_sample)
rT_est_termo=np.polyval(termo_fit, log_k_sample) # stima del reciproco della T dato il log (K_sample)
# sulla base del fit appena fatto (i parametri a e b
# del fit sono nella variabile termo_fit)
T_est_termo=1/rT_est_termo # passiamo dal "reciproco della T" a T
print("Temperatura stimata: %5.2f\n" % T_est_termo)
La differenza tra le due stime di temperatura è circa 40K. Non è una grandissima differenza, ma è pur sempre ben maggiore dell'errore che ci aspettiamo sulla base di questo tipo di misura, come abbiamo visto nella prima parte del notebook.
Se avessimo avuto tutti i dati (anche quelli di più bassa temperatura), avremmo ottenuto una stima molto vicina a quella che abbiamo dato usando il fit corretto:
temperatura(0.010,0.162)
Morale: studiare termodinamica!