\chapter{Equazione del calore}

\section{Impostazione}
Supponiamo che ad un corpo di volume $V$ venga fornita {\em in un determinato istante} una quantit\`a di calore $Q$. Nel volume elementare $dV$ ci sar\`a  una variazione di quantit\`a  di calore $dQ=C_v dm dT$, essendo $C_v$ la capacit\`a termica a volume costante, $dm=\rho dV$ la massa del volume $dV$ con densit\`a  $\rho$ e $dT$ la variazione di temperatura. Dunque, nell'unit\`a  di tempo, sar\`a:
\begin{equation}
(dQ)_v = C_v \rho dV \frac{\partial T}{\partial t} dt
\label{eq:var_calore}
\end{equation}
Il flusso di calore che passa attraverso la superficie $dS$ del nostro volume elementare $dV$ nell'unit\`a  di tempo è:
\begin{equation}
(dQ)_s = k dS dt \frac{\partial T}{\partial \mathbf n}
\label{eq:flusso_elem}
\end{equation}
essendo $k$ la conducibilit\`a termica (dimensioni $[k]=W m^{-1} K^{-1}$), e 
$$\frac{\partial T}{\partial \mathbf n}=\nabla T\cdot \mathbf n$$
la derivata dela temperatura nella direzione normale alla superficie considerata.\par
Integrando su tutta la superficie del volume $V$, si ottiene, usando la formula di Ostrogradskij:
\begin{equation}
dQ=\int_{S} k dt \nabla T \cdot n dS=\int_{V} k dt \nabla^2 T dV
\label{eq:flusso_tot}
\end{equation}
e dunque, uguagliando alla variazione di calore (\ref{eq:var_calore}):
\begin{equation}
\int_{V} \left( \rho C_v \frac{\partial T}{\partial t} - k \nabla^2 T \right) dV =0
\end {equation}
Siccome questa equazione deve essere valida per qualsiasi volume, possiamo eliminare l'integrazione e scrivere l'equazione della trasmissione del calore:

\begin{equation}
\boxed{\frac{\partial T}{\partial t} - a^2 \nabla^2 T = 0}
\label{eq:trasm_calore}
\end{equation}
avendo definito la diffusivit\`a termica $a^2=\frac{k}{\rho C_v}$, con dimensioni $a^2=m^2 s^{-1}$.

\subsection{Geoterme e lunghezza di diffusione}
In presenza di una sorgente distribuita di calore nel volume considerato, se indichiamo con $A$ la densità delle sorgenti, l'equazione di diffusione del calore diventa:
\begin{equation}
k \nabla^2 T+A=\rho C_v \frac{\partial T}{\partial t}
\end{equation}
che, in una situazione di stato stazionario -- ovvero quando la produzione di calore e la sua dispersione sono in equilibrio e non c'è variazione di temperatura nel tempo -- diventa, semplicemente:
\begin{equation}
\nabla^2 T= \frac{A}{k}
\label{eq:geoterme}
\end{equation}
L'eq.(\ref{eq:geoterme}), di permette di stabilire la variazione della temperatura con la profondità. La dipendenza funzionale della temperatura con la profondità prende il nome di {\em geoterma}. Nel caso di assenza di sorgenti, la temperatura varierà linearmente con la profondità. Infatti, per $A=0$ e considerando solo la direzione z:
\begin{displaymath}
\nabla^2 T=0 \Rightarrow \frac{\partial T}{\partial z}= cost = C \Rightarrow T=Cz+T(z=0)
\end{displaymath}
Un valore tipico per la costante C (la geoterma) è di $20 K km^{-1}$, che, combinato con il valore medio della conduttività $k=3.0 W m^{-1} K^{-1}$, fornisce il flusso di calore per unità di superficie terrestre: $0.06 W m^{-2}$. Con questi valori, per una profondità di 60 km, si raggiungerebbe una temperatura di $20 \times 60 + 300=1500 K$, che è superiore alla temperatura di fusione della maggior parte delle rocce. D'altra parte, da misure di propagazione di onde di taglio nel mantello, risulta che questo si comporta (su scale di tempi brevi) come un solido.\par
In presenza di sorgenti distribuite nel volume -- come per es. in presenza di materiali radioattivi che forniscono calore per decadimento -- sarà:
\begin{displaymath}
\frac{\partial T}{\partial z}=\frac{A}{k}x + C \Rightarrow T=\frac{A}{2k}x^2 + C x +D
\end{displaymath}
che, a profondià maggiori, peggiorerebbe la discrepanza con i dati sismici. Bisogna, però, tenere conto di due fattori:
\begin{itemize} 
\item La temperatura di fusione aumenta con la pressione, e di conseguenza con la profondità
\item Oltre una certa profondità, il gradiente di temperatura non è più determinato dal raffreddamento conduttivo ma dalla compressione adiabatica. Il gradiente dovuto a questo processo, in cui la temperatura aumenta per sola compressione, senza scambio con l'esterno, è molto minore di quello conduttivo.
\end{itemize}
\subsubsection{Lunghezza di diffusione}
L'analisi dimensionale del coefficiente di diffusività, $a^2=\frac{k}{\rho C_v}$ suggerisce una lunghezza di scala $L$ -- la {\em lunghezza di diffusione} -- tale che una variazione di temperatura all'istante $t_0$, all'istante $t$, si sia propagata ad una distanza L nel mezzo la cui diffusività è $a^2$.\par
Per verificarlo, consideriamo variazioni periodiche di temperatura alla superficie:
\begin{displaymath}
T(0, t)=T_{max}e^{-i\omega t}
\end{displaymath}
 Ad una profondità $z$, la temperatura avrà l'andamento:
\begin{equation}
T(z, t)=T_{max}e^{i(kz-\omega t)}
\label{eq:temp_perio}
\end{equation}
Sostituendo nell'equazione (\ref{eq:trasm_calore}), otteniamo la relazione {\em di dispersione}
\begin{displaymath}
a^2k^2=i\omega
\end{displaymath}
da cui dicaviamo il "numero d'onda" 
\begin{equation}
k=\frac{\sqrt{i\omega}}{a}
\end{equation}
che, sostituito nella (\ref{eq:temp_perio}) dà:
\begin{equation}
T(z, t)=T_{max}e^{\frac{i^{\frac{3}{2}}\omega^{\frac{1}{2}}}{a}z-i\omega t}
\end{equation}
Ricordando che $i^{\frac{3}{2}}=\frac{i-1}{\sqrt{2}}$, raccogliendo i termini reali e immaginari all'esponente:
\begin{equation}
T(z, t)=T_{max}e^{-\frac{z}{\delta}}e^{i(\frac{z}{\delta} -\omega t)}
\end{equation}
dove abbiamo definito la lunghezza di penetrazione (o {\em skin depth}) 
\begin{equation}
\delta=\sqrt{\frac{2}{\omega}}a
\end{equation}
come la profondità alla quale l'ampiezza delle oscillazioni si è ridotta ad un fattore 1/e di quella iniziale.\par
La skin depth, introduce anche un ritardo di fase $\Delta \phi$ delle oscillazioni visibili alla profondità $z$ rispetto a quelle in superficie:
\begin{equation}
\Delta \phi = \frac{z}{\delta}=\left(\frac{\omega}{2a^2}\right)^{\frac{1}{2}}
\end{equation}
Ad esempio, per una profondità di 1 m d'acqua ($a^2= 1.4 \cdot 10^{-7} m^2 s^{-1}$), il ciclo giorno-notte
($\omega = \frac{2 \pi}{24 h\times 3600 s/h}=7.3\times 10^{-5} s^{-1}$), è spostato in avanti di 16 secondi.


\section{Soluzione nel caso monodimensionale}
Come nel caso dell'equazione delle onde possiamo cercare la soluzione dell'equazione monodimensionale della trasmissione del calore, con la fattorizzazione della soluzione in una parte $\tau(t)$, dipendente dal tempo ed una $X(x)$, dipendente dallo spazio. L'equazione (\ref{eq:trasm_calore}):
\begin{equation}
T_t-a^2T_{xx}=0
\end{equation}
con la condizione iniziale: 
\begin{equation}
T(x,0)=f(x)
\label{eq:cal_cond}
\end{equation}
ponendo $T(x, t)= X(x)\tau(t)$, diventa:
\begin{equation}
\dot\tau(t)X(x)-a^2 \tau(t)X''(x) = 0
\end{equation}
Come nel caso dell'eq. delle onde, con le stesse considerazioni per la parte spaziale, possiamo dire:
\begin{equation}
\frac{\dot\tau(t)}{a^2\tau(t)}=\frac{X''(x)}{X(x)} =-\lambda^2
\label{eq:calor_separata}
\end{equation}
che, per la parte temporale d\`a, immediatamente:
\begin{equation}
\tau(t)=C e^{-a^2 \lambda^2 t}
\end{equation}
 In totale, la soluzione si pu\`o scrivere come:
 \begin{equation}
 T(x,t)=\int_{-\infty}^{\infty} C(\lambda) e^{-a^2\lambda^2 t} e^{i\lambda x} d\lambda
 \label{eq:calor_sol_compl}
\end{equation}
La condizione iniziale (\ref{eq:cal_cond}), pu\`o essere scritta esprimendo la f(x) in funzione della sua trasformata di Fourier (\ref{eq:trasf_fourier}):
 \begin{align}
 \nonumber
T(x,0)=f(x)=&\frac{1}{2\pi}\int_{-\infty}^{\infty} F(\lambda)e^{i\lambda x} d\lambda=\\
\nonumber
&=\frac{1}{2\pi}\int_{-\infty}^{\infty} \left[\int_{-\infty}^{\infty}f(x')e^{-i\lambda x'}dx'\right] e^{i\lambda x} d\lambda=\\
\nonumber
&= \int_{-\infty}^{\infty} C(\lambda) e^{i\lambda x} d\lambda
\end{align}
che è verificata per 
\begin{equation}
C(\lambda)=  \frac{1}{2\pi}\int_{-\infty}^{\infty}f(x')e^{-i\lambda x'} dx'
\label{eq:calor_costX}
\end{equation}
Sostituendo nella soluzione (\ref{eq:calor_sol_compl}) otteniamo:
\begin{eqnarray}
\nonumber
\frac{1}{2\pi}\int_{-\infty}^{\infty}f(x')dx'
\left[\int_{-\infty}^{\infty}e^{-a^2\lambda^2 t}e^{i\lambda (x-x')}d\lambda\right]=\\
=\frac{1}{\pi}\int_{-\infty}^{\infty}f(x')dx\int_{0}^{\infty}e^{-a^2\lambda^2 t}cos(\lambda(x-x'))d\lambda
\label{eq:sol_interm}
\end{eqnarray}
avendo sfruttato la simmetria dell'integrando rispetto a $\lambda$ e la definizione di coseno in termini di esponenziali complessi. L'integrale (\ref{eq:sol_interm}) \`e del tipo (Smirnov [81])
$$
\int_{0}^{\infty}e^{-\alpha^2 x^2}cos(\beta x) dx = \frac{\sqrt{\pi}}{2\alpha}e^{-\frac{\beta^2}{4\alpha^2}}
$$
Dunque, la soluzione completa dell'equazione di trasmissione del calore sar\`a:
\begin{equation}
T(x, t)=\int_{-\infty}^{\infty}f(x')\frac{1}{2a\sqrt{\pi t}}e^{-\frac{(x-x')^2}{4a^2t}}dx'
\label{eq:calor_sol_compl1}
\end{equation}
In particolare, se $f(x')=\delta(x')$ (il calore inizialmente è somministrato in un punto a $x'$), la soluzione \`e l'integrando che \`e la distribuzione di temperatura risultante dalla somministrazione istantanea all'istante $t=0$ di una quantit\`a di calore $\rho C_v$ nel punto x=x'.


\section{Soluzione per spazio semi-infinito}
\label{sec:calor_semi_inf}
La soluzione dell'equazione si pu\`o trovare pi\`u facilmente considerando uno spazio semi-infinito, come per es. il caso della propagazione del calore dalla superficie terrestre verso l'interno (o voceversa) quando si consideri infinito il raggio terrestre rispetto alle distanze di interesse ([Turcotte 4.15]).\\
Riscriviamo l'equazione in termini della profondit\`a:
$$
\frac{\partial T}{\partial t}=\frac{k}{\rho C_v } \frac{\partial^2 T}{\partial z^2}
$$
con le condizioni iniziali:
\begin{equation}
\left\{ 
\begin{array}{l}
T(z, t=0)=T_0\\
T(z=0, t)=T_1\\
T(z=\infty, t)=T_0
\end{array}
\right.
\label{eq:semiinf_cond}
\end{equation}
che indicano che all'istante t=0, il semispazio, inizialmente  \`e a temperatura $T_0$, per t>0, la superficie viene mantenuta a temperatura $T_1$. Il calore si propagher\`a verso l'interno se $T_1>T_0$ (viceversa il semispazio si raffredder\`a se $T_1<T_0$) e a profondit\`a infinita, la temperatura sar\`a quella iniziale.
Se facciamo il cambiamento di variabile:
\begin{equation}
\theta=\frac{T(z,t)-T_0}{T_1-T_0}
\label{var_theta}
\end{equation}
l'Equazione diventa:
\begin{equation}
\frac{\partial \theta}{\partial t}=\alpha\frac{\partial^2 \theta}{\partial x^2}
\label{eq:semii_intheta}
\end{equation}
le condizioni(\ref{eq:semiinf_cond}), saranno:
\begin{equation}
\left\{ 
\begin{array}{l}
\theta(z, 0) = 0\\
\theta(0,t) = 1\\
\theta(\infty, t)=0
\end{array}
\right.
\label{eq:semiinf_cond1}
\end{equation}
Effettuando, l'ulteriore cambiamento di variabili :
\begin{equation}
\eta(z,t)=\frac{z}{2\sqrt{\alpha t}}
\label{var:eta}
\end{equation}
risulta:
\begin{eqnarray}
\nonumber
\frac{\partial \theta}{\partial t}= \frac{\partial \theta}{\partial \eta}\frac{\partial \eta}{\partial t}=\\
\nonumber
= \frac{\partial \theta}{\partial \eta}\left(-\frac{z}{2\sqrt{\alpha t}}\frac{1}{2t}\right)=
 -\frac{\eta}{2t} \frac{\partial \theta}{\partial \eta}
\end{eqnarray}
\begin{eqnarray}
\nonumber
\frac{\partial \theta}{\partial z}= \frac{1}{2\sqrt{\alpha t}}\frac{\partial \theta}{\partial \eta}\\
\frac{\partial^2 \theta}{\partial z^2}=\frac{1}{2\sqrt{\alpha t}}\frac{\partial^2 \theta}{\partial \eta^2}\frac{\partial \eta}{\partial z}=\frac{1}{4\alpha t}\frac{\partial^2 \theta}{\partial \eta^2}
\end{eqnarray}
Sostituendo nella (\ref{eq:semii_intheta}), otteniamo:
\begin{equation}
-\eta\frac{\partial{\theta}}{\partial \eta}=\frac{1}{2}\frac{\partial^2 \theta}{\partial \eta^2}
\label{eq:intheta_e_eta}
\end{equation}
che \`e nella sola variabile $\eta$ e pu\`o essere facilmente integrata.\\
Infatti, se $\phi=\frac{\partial \theta}{\partial \eta}$, allora: 
$$-\eta \phi=\frac{1}{2}\frac{\partial \phi}{\partial \eta}\Rightarrow \phi=c_1e^{(-\eta^2)}$$
Dunque:
$$\frac{\partial\theta}{\partial \eta}= c_1 e^{-\eta^2}\Rightarrow \theta=\int_{0}^{\eta}c_1 e^{-\eta'^2}d\eta'+c_2$$
 le costanti $c_1$ $c_2$ si ricavano dalle condizioni iniziali $\theta(\infty)=0$ e $\theta(0)=1$.  
 Dalla prima condizione viene:
 \[
 c_1\int_{0}^{\infty} e^{-\eta'^2}d\eta=-c_2 = c_1\sqrt{\pi}/2
 \]
 e dalla seconda:
\[
c_2 =1
\]
 Dunque risulta:
\begin{equation}
\theta = 1-\frac{2}{\sqrt{\pi}}\int_{0}^{\eta} e^{\eta'^2}d\eta'=1-Erf(\eta)=Erfc(\eta)
\end{equation}
\begin{figure}[t]
\begin{center}
\includegraphics[width=0.7\textwidth]{Fig_Erf.eps}
\caption{Andamento funzionale delle funzioni Erf($\eta$) e Erfc($\eta$) per $\eta$ positivi}
\end{center}
\label{Fig:Erf}
\end{figure}
Ritornando alle variabili iniziali:
\begin{eqnarray}
\nonumber
\frac{T-T_0}{T_1-T_0}=1-Erf(\eta) \Rightarrow (T_1-T_0)-(T_1-T_0)Erf(\eta)+T_0=T(z, t)\\
\boxed{T(z,t)=T_1\left[1-Erf\left(\frac{z}{2\sqrt{\alpha t}}\right)\right]+ T_0Erf\left(\frac{z}{2\sqrt{\alpha t}}\right)}
\label{eq:sol_calor_bidim}
\end{eqnarray}
Per $t\rightarrow\infty$ o $z=0$, vale $T(z, t=\infty)= T(z=0,t)=T_1$ (che significa che per tempi sufficientemente lunghi tutto il sistema si sarà portato alla temperatura superficiale), mentre, per tempi finiti e grandi profondità $T(z\rightarrow\infty, t)=T_0$.
\subsection{Flusso di calore}
Se conosciamo la temperatura in superficie, possiamo calcolare il flusso di calore dall'interno della terra, misurando il gradiente.\par
Il flusso è:
\begin{align*}
\Phi_u=&-k \left.\frac{\partial T}{\partial z}\right|_{z=0}=\\
& = -k\left(T_{sup} - T_0 \right)\frac{\partial}{\partial z}\left[ 1-erf\left(\frac{z}{2\sqrt{\alpha t}}\right)\right]_{z=0}=\\
& = k \left. \left( T_{sup} - T_0 \right) \frac{1}{2\sqrt{\alpha t}}\frac{2}{\sqrt{\pi}} e^{-\eta^2}\right|_{\eta=0} 
\end{align*}
Dunque, per $t>0$ sarà:
\begin{equation}
\Phi_u= k \frac{\left(T_{sup} - T_0\right)}{\sqrt{\pi\alpha t}}
\label{eq:therm_flux}
\end{equation}
Da cui si vede che il flusso di calore decresce con l'inverso della radice dell'età della litosfera, e dunque si può determinare l'età della litosfera misurando il flusso.\par
\subsubsection{Stima dell'età della Terra}
William Thomson (in seguito noto come Lord Kelvin), stimò l'età della Terra assumendo che la perdita di calore avvenisse per sola conduzione. Supponiamo che $T_{sup}=0$ e $T_0\sim 2\cdot 10^3$ $ ^{o}C$, e sia $\alpha=10^{-6}m^2 s^{-1}$. Scrivendo la (\ref{eq:therm_flux}):
\begin{displaymath}
k \left.\frac{\partial T}{\partial z}\right|_{z=0}= \frac{k T_0}{\sqrt{\pi\alpha t}}
\end{displaymath} 
da cui:
\begin{displaymath}
t=\frac{T_0^2}{\pi \alpha\left(\frac{\partial T}{\partial z}\right)_{z=0}}
\end{displaymath}
Assumendo un gradiente di temperatura (misurato) di $3\cdot 10^{-3}$ $^o C m^{-1}$, si otterrebbe un valore per l'età della Terra di $\sim 4.4 \cdot 10^7$ anni che è più di un ordine di grandezza in meno di quanto ci si aspettava già all'epoca in base a numerose altre prove. Da questo si dedusse che doveva esserci qualche altra sorgente di calore oltre al nucleo. Successivamente, questa fonte, venne attribuita alla presenza di sostanze radioattive distribuite all'interno della crosta, che perdono energia per decadimento radioattivo e forniscono il flusso mancante per stimare una giusta età. In realtà, anche questa sorgente non basta a spiegare l'età della terra, giacché la produzione di calore per decadimento radioattivo, si stima sia circa la metà dell'odierno flusso di calore. La soluzione, in realtà, fu data circa 30 anni dopo la stima di Thomson da John Perry, che, molto semplicemente, mostrò che la variazione della conduttività con la temperatura avrebbe cambiato di molto il risultato. In particolare, mostrò che il trasporto di calore per convezione avrebbe alterato grandemente il risultato di Thomson.
\subsubsection{Perdita di calore per decadimento}
Il meccanismo del decadimento radioattivo deriva dalla circostanza che alcuni nuclei atomici hanno una probabilità finita di perdere neutroni o particelle $\alpha$ (nuclei di elio), trasformandosi in nuclei più leggeri e cedendo all'ambiente la differenza di energia tra lo stato legato iniziale e quello finale.\par
Se prendiamo un certo numero di questi atomi, e verifichiamo, nell'unità di tempo quanti di questi sono decaduti, questa frazione ci darà la probaliltà di decadimento $P$ del singolo atomo. In altri termini, preso un atomo, in un tempo $\tau$, questi avrà una probabilità $P=1/\tau$ di decadere. Questo ci fornisce la legge di decadimento. Infatti, se $N_0$ è il numero di atomi non decaduti iniziali e il numero di atomi al tempo $t$ è $N(t)$, allora il numero di atomi non decaduti al tempo $t+dt$ sarà:
$$ N(t+dt)=N(t)+\frac{dN}{dt}dt=N(t)\left(1-\frac{dt}{\tau}\right)$$
avendo supposto che in un tempo $dt$ siano decaduti $P\cdot N(t)dt$ atomi. Dunque:
\begin{displaymath}
\frac{dN}{N}=-\frac{dt}{\tau}\Rightarrow N(t)=N_0e^{-\frac{t}{\tau}}
\end{displaymath}
Dunque possiamo definire un tempo -- il tempo di dimezamento -- dopo il quale gli atomi della specie iniziale si sono ridotti alla metà, dunque quando $dN=1/2N$ e quindi $d(ln N)= 1/2$:
$$ N(t)=N_0 e^{-\frac{t_{1/2}}{\tau}}=\frac{1}{2}N_0\Rightarrow t_{1/2}=\tau ln(2)$$
ed esprimere i tempi in termini del tempo di dimezzamento.\par
La concentrazione di atomi di una determinata specie radioattiva $x$ ad un epoca t, sarà dunque funzione del tempo di dimezzamento caratteristico della specie:
\begin{equation}
C^x(t)=C_0^xe^{-\frac{t ln2}{\tau_{1/2}}}
\end{equation}
La quantità di calore $H$ rilasciata dal complesso delle sostanze radioattive presenti nella crosta terrestre, sarà data dalla somma delle quantità di calore rilasciate da ciascuna specie, moltiplicate per la relativa concentrazione. I principali componenti che contribuiscono al flusso di calore nella crosta terrestre sono due isotopi dell'Uranio ($^{238}U$ e $^{235}U$) il Torio ($^{232}Th$) e un isotopo del potassio ($^{40}K$). Nella tabella che segue, diamo le concentrazioni, le quantità di calore rilasciate e i tempi di dimezzamento degli isotopi citati:\\
\begin{center}
\begin{tabular}{|c|c|c|c|} 
\hline
 {\bf Elemento} & {\bf H}  & $\tau_{\frac{1}{2}}$& {\bf Concentrazione}\\
\hline
$^{238}U$ &  $9.37 \cdot 10^{-5}W kg^{-1}$ & $4.4 \cdot 10^8$ anni & 99.27\% U\\
\hline
$^{235}U$ & $1.36\cdot 10^{-7}W kg^{-1}$ &  $7.04 \cdot 10^8$ anni & 0.73\% U\\
\hline
$^{232}Th$ & $2.69 \cdot 10^{-5}W kg^{-1}$ & $1.4 \cdot 10^{10}$anni & 100 \% Th\\
\hline
$^{40}K$ & $2.79 \cdot 10^{-5}W kg^{-1}$ &$1.25 \cdot 10^9$ anni & 0.0128 \% K \\
\hline 
\end{tabular}
\end{center}
Con questi dati, la produzione di calore dovuta alle sostanze radioattive, in funzione del tempo si può scrivere come:
\begin{align}
\nonumber
H =& 0.9927 C_0^U H_{^{238}U} e^{-\frac{t ln2}{\tau_{1/2}^{^{238}U}}}+ 0.0073C_0^U H_{^{235}U} e^{-\frac{t ln2}{\tau_{1/2}^{^{235}U}}}\\
&+ C_0^{Th} H_{Th} e^{-\frac{t ln2}{\tau_{1/2}^{Th}}}+ 1.28\cdot10^{-4}C_0^K H_{^{40}K} e^{-\frac{t ln2}{\tau_{1/2}^{^{40}K}}}
\label{eq:decay}
\end{align}
L'andamento nel tempo della produzione di calore dovuta al decadimento radioattivo basato sulla (\ref{eq:decay}) è mostrato in fig. \ref{fig:decay}, da cui si evince che, a causa della vita media più corta, nel passato $^{235}U$ e $ ^{40}K$, erano dominanti, mentre ora predominano $^{238}U$ e $^{232}Th$.\par
\begin{figure}[ht]
\begin{center}
	\includegraphics[width=0.7\textwidth]{decay.pdf}
	\caption{Andamento nel tempo della produzione di calore dovuto al decadimento degli isotopi radioattivi e contributo di ciascun isotopo}
	\end{center}
	\label{fig:decay}
\end{figure}
Consideriamo ora una distribuzione uniforme $A=\rho H_0$ di sorgenti radiogeniche all'interno della crosta continentale ed uno spessore medio dei continenti di 35 km. Allora sarà:
\begin{displaymath}
\frac{\partial^2 T}{\partial z^2}= \frac{\rho H_0}{k}
\end{displaymath}
con le condisioni iniziali:
\begin{displaymath}
 \left\{
\begin{array}{rl}
T(z=0)=&T_0\\
-k\frac{\partial T}{\partial z}=& -\Phi_u
\end{array}
\right.
\end{displaymath}
e dunque:
\begin{displaymath}
T(z)=T_0+\frac{\Phi_u z}{k}-\frac{\rho H_0}{2 k} z^2
\end{displaymath}
Se consideriamo la composizione della crosta continentale come quella del granito, possiamo porre:
$H_0^{gr}=1.6\cdot 10^{-10}\frac{W}{kg}$ e $\rho^{gr}=2700 \frac{kg}{m^3}$, dunque:
\begin{displaymath}
\left. k\frac{\partial T}{\partial z}\right|_{z=D}=\Phi_u - \rho H_0 D
\end{displaymath}
Siccome a $z=D$ la derivata deve essere nulla:
\begin{displaymath}
	\Phi_u = \rho H_0 D = 2.7\cdot 10^3\cdot3.5\cdot 10^4\cdot 9.6\cdot 10^{-10}=91\, mW/m^2
\end{displaymath}
che è maggiore del flusso realmente osservato di $65\, mW/m^2$.


\subsection{Spessore del fondo oceanico}
\begin{figure}[t]
\begin{center}
\includegraphics[width=0.7\textwidth]{crosta_oceanica.eps}
\caption{Modello per la descrizione dello spessore della crosta oceanica in funzione della distanza dal ridge}
\end{center}
\label{Fig:Crosta}
\end{figure}
Il modello delineato in Fig. \ref{Fig:Crosta}, rappresenta la metà di un fondo oceanico in espansione. Le assunzioni del modello sono che:
\begin{itemize}
\item Le placche litosferiche sono rigide e si allontanano dalla frattura medio-oceanica (il ridge) con velocità costante $v$.
\item Il vuoto viene riempito da materiale caldo a bassa viscosità dall'astenosfera
\item la gnerazione di calore interna è molto inferiore agli altri termini dell'equazione del calore e può essere trascurata
\item c'è una singolarità a $x=z=0$
\end{itemize}
Si tratta, di un problema bidimensionale senza sorgenti e quindi l'equazione sarà:
\begin{equation}
\frac{\partial ^2 T}{\partial x^2}+\frac{\partial ^2 T}{\partial z^2}=\frac{v}{k}\frac{\partial T}{\partial x}
\label{eq:expa_ridge}
\end{equation}
dove al secondo membro c'è il termine advettivo, ovvero il termine di trasporto di calore dovuto al movimento della placca. Si può mostrare che il termine di diffusione laterale si può trascurare per $x>>0$. Considerato che possiamo scrivere: 
\begin{displaymath}
\frac{\partial T}{\partial x}=\frac{\partial T}{\partial t}\frac{\partial x}{\partial t}
\end{displaymath}
l'eq. (\ref{eq:expa_ridge}) diventa la solita equazione del calore.\par
Per stabilire lo spessore della crosta oceanica dobbiamo definire quale sarà la sua base, ovvero dove la crosta da solida passerà ad essere liquida. Consideriamo che questo passaggio avvenga ad una temperatura che è una frazione abbastanza grande della temperatura del mantello $T_m$, diciamo il 90\%. La soluzione del problema bidimensionale l'abbiamo già vista in Eq.(\ref{eq:sol_calor_bidim}) e dunque possiamo ricavare la profondità:
\begin{eqnarray}
\nonumber
\frac{T(z,t) - T_m}{T_{sup}-T_m}=&\left(1-Erf\left(\frac{z}{2\sqrt{\alpha t}}\right)\right)\\
\end{eqnarray}
Se imponiamo che $T(z,t)=0.9 T_m$ e trascuriamo $T_{sup}$ rispetto a $T_m$ otteniamo:
\begin{equation}
0.9 T_m= T_m Erf\left(\frac{z}{2\sqrt{\alpha t}}\right)
\end{equation}
La profondità cercata è, dunque quella per cui $Erf(\eta)=0.9$. I valori della funzione inversa della $Erf(\eta)$, si posssono trovare tabellati nei testi di statistica (esiste anche una funzione di Matlab), e si trova: $\eta|_{T=0.9T_m}=Erf^{-1}(0.9)=1.1631$. Di conseguenza: 
\begin{displaymath}
z|_{t=0.9T_m}=2\eta_T \sqrt{\alpha t}=2.326 \sqrt{\alpha t}
\end{displaymath}
Se consideriamo un valore di $\alpha=8\cdot 10^{-7} m^2 s^{-1}$ e che in un anno ci sono $3\cdot 10^7$ secondi, otteniamo un valore di profondità di $z(km)\sim 11.4 \sqrt{t (Ma)}$.\par
\subsubsection{Profondità del fondale oceanico}
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.9\textwidth]{depth.eps}
\caption{Modello per la descrizione della profondità oceanica in funzione del tempo. L'asse delle profondità non è in scala}
\end{center}
\label{Fig:Eq_hydro}
\end{figure}
Stabilito che lo spessore della crosta oceanica aumenta con l'età e, di conseguenza, con la distanza dal ridge, siccome la crosta oceanica è in equilibrio idrostatico (galleggia) con il mantello, ne deriva che la profondità del fondo oceanico aumenta con la distanza dal ridge. Per dare un'espressione più quantitativa, consideriamo anche l'andamento della densità della crosta in funzione della profondità $\rho(z,t)$. Dette $\rho_w=1.0 \cdot 10^3 kg \, m^{-3}$ la densità dell'acqua e $\rho_m=3.3  \cdot 10^3 kg \, m^{-3}$ quella del mantello, e fissata a 0 la profondità del ridge, l'espressione dell'equilibrio idrostatico di una sezione di litosfera spessa $z_L$ a profondità $w$ sarà(vedi Fig. \ref{Fig:Eq_hydro}):
\begin{equation}
(w+z_L)\rho_m=w\rho_w+\int_0^{z_L}\rho(z, t) \,dz
\label{eq:eq_hydro}
\end{equation}
Raccogliendo i termini in $w$, la (\ref{eq:eq_hydro}) si può scrivere:
\begin{equation}
(\rho_w-\rho_m)w+\int_0^{z_L}(\rho(z,t)-\rho_m)dz = 0
\label{eq:eq_hydro1}
\end{equation}
da cui si vede che, visto che il primo termine è negativo, il secondo deve essere positivo, ovvero $\rho(z,t) >  \rho_m$. I cambiamenti di pressione e temperatura con la profondità, producono dei cambiamenti di volume e quindi di densità:
\begin{equation}
dV=\left(\frac{dV}{dT}\right)_P + \left(\frac{dV}{dP}\right)_T
\end{equation}
Comunque, fissata una posizione ad una certa distanza dal ridge, non ci sono variazioni di pressione, e dunque rimane solo il primo termine. D'altra parte il primo termine è legato al coefficiente di espansione termica 
$$\alpha = \frac{1}{V}\left(\frac{dV}{dT}\right)_P$$
e, visto che $\rho=m V^{-1}$, implica
$$\frac{d\rho}{\rho}=-\frac{dV}{V}$$
risulta:
\begin{equation}
d\rho=-\alpha \rho dT
\label{eq:therm_exp}
\end{equation}
e, quindi, applicandola al nostro caso:
\begin{displaymath}
\rho(z,t) - \rho_m=-\alpha\rho_m [T - T_m]
\end{displaymath} 
che, sostituita nella (\ref{eq:eq_hydro1}), ci da:
\begin{equation}
(\rho_w - \rho_m)w - \rho_m \alpha\int_0^{z_L} [T_m - T(z,t)] dz=0
\end{equation}
D'altra parte, dalla soluzione dell'equazione del calore, sappiamo che:
$$ \frac{T-T_m}{T_s - T_m}=1-Erf(\eta) $$
e dunque 
\begin{equation}
(\rho_m - \rho_w)w = \alpha \rho_m [T_m - T_s]\int_0^{z_L}Erfc\left(\frac{z}{2\sqrt{at}}\right) dz
\end{equation}
Possiamo, con buona approssimazione, estendere l'integrale fino all'infinito ottenendo:
\begin{equation}
(\rho_m - \rho_w)w = \alpha \rho_m [T_m - T_s] 2\sqrt{\frac{at}{\pi}}
\end{equation}
da cui  si ricava l'espressione della profondità oceanica in funzione dell'età:
\begin{equation}
\boxed {w(t)=\sqrt{\frac{at}{\pi}}\frac{2\alpha \rho_m}{\rho_m - \rho_w}(T_m - T_s) }
\label{eq:Oc_depth}
\end{equation}
che ci dice che la profondità degli oceani cresce con la radice dell'età della litosfera, ovvero con la radice della distanza dal ridge.
Diamo ora una tabella con i numeri per stimare la profondità:
\begin{center}
\begin{tabular}{|c|c|c|} 
\hline
Parametro & definizione & valore \\
\hline
$T_s$ & Temperatura superficiale & $300 \, K$\\
\hline
$T_m$ & Temperatura del mantello & $1600 \, K$ \\
\hline
$a$ & Diffusività termica & $8 \cdot 10^{-7} m^2 s^{-1}$ \\
\hline
$\alpha$ & Coefficiente di espansione termica & $3.1 \cdot 10^{-5}\, K^{-1}$\\
\hline
$\rho_w$ & Densità dell'acqua marina & $1025\, kg\, m^{-3}$ \\
\hline
$\rho_m$ & Densità del mantello & $3300\, kg\, m^{-3}$\\
\hline
\end{tabular}
\end{center}
Con questi valori, tenendo conto che $1\, anno\sim\pi\cdot 10^7 s$, si ottiene che la profondità degli oceani è di circa $d(m)=350 \sqrt{t(Ma)}$ al di sotto della profondità dell'asse del ridge. 
