Consider the following model for fitting: \[P(t)=ae^{\lambda (t-2000)}\] where a and λ are the paramters we seek to fit. Apply a simple linear regression on the previous graph and we get:
| RES | a | λ |
|---|---|---|
| Hydro | 304.8085830 | 0.0209350 |
| Wind | 19.0710527 | 0.2270686 |
| Solar | 0.5769026 | 0.4231368 |
| Geothermal | 8.4545894 | 0.0607396 |
We fit the model into the original graph for comparison. λ is the growth rate of each renewable source.
A higher λ implies a higher growth rate for the RES. It is mostly affected by the learning curve of the RES. For example, costs in deployment has dropped significantly for wind and solar in the previous years, thus yielding a much higher λ than hydro and geothermal.
For constant radiogenic sources we have \[q(z)=-(q_{s}-(q_{s}-q_{b}) \dfrac{z}{h});z<h\] \[q(z)=-q_{b};z>h\]
from Fourier’s Law we have \[\dfrac{\partial T}{\partial z}=-\dfrac{q(z)}{\lambda}\]
thus the temperature profile is \[T(z)=T_{s}+\dfrac{q_{s}z-(q_{s}-q_{b}) \dfrac{z^2}{2h}}{\lambda};z<h\] \[T(z)=T_{s}+\dfrac{(q_{s}+q_{b}) \dfrac{h}{2}+q_b(z-h)}{\lambda};z>h\]
For exponentially decaying radiogenic sources we have \[q(z)=-(q_{b}+(q_{s}-q_{b})e^{-\dfrac{z}{h}})\]
Similarly we can deduce \[T(z)=T_{s}+\dfrac{q_{b}z+(q_{s}-q_{b})h(1-e^{-\dfrac{z}{h}})}{\lambda}\]
Below we plot the heat fluxes for constant cases (lines) and exponential decaying cases (dots) for different \(q_s\).
Below we plot the geotherms for constant cases (lines) and exponential decaying cases (dots) for different \(q_s\).
For constant radiogenic sources, we have the depth for certain level of temperature \[z=\dfrac{h}{1-\dfrac{q_b}{q_s}}(1-\sqrt{1-\dfrac{2\lambda(1-\dfrac{q_b}{q_s})(T_z-T_s)}{hq_s}});z<h\] \[z=h+\dfrac{1}{q_b}(\lambda(T_z-T_s)-\dfrac{(q_s+q_b)h}{2});z>h\]
The depth for certain level of temperature for the exponentially decaying case cannot be solved explicitly. Instead, we solve for the relationship between surface heat flux and the depth reaching that temperature: \[(\dfrac {\partial z}{\partial q_s})_{T}=-\dfrac{(\dfrac{\partial T}{\partial q_s})_{z}}{(\dfrac{\partial T}{\partial z})_{q_s}}\]
for exponentially decaying radiogenic sources this is \[(\dfrac {\partial z}{\partial q_s})_{T}=\dfrac{h(1-e^{-\dfrac{z}{h}})}{q(z;q_s)}\]
With an initial knowledge of any given point that satisfies \(T(z;q_s)=T_z\), we can solve the ODE with numerical methods. We use RK-4 here.
RK<-function(z,q_s,dq){
k_1=delta(z,q_s)
k_2=delta(z+.5*k_1*dq,q_s)
k_3=delta(z+.5*k_2*dq,q_s)
k_4=delta(z+k_3*dq,q_s)
z_new=z+(k_1+2*k_2+2*k_3+k_4)*dq/6
return(z_new)
}The resulting \(z-q_s\) plot is
The exponential decaying model assumes a greater amount of radiogenic source than the constant source model. Therefore, the heat flux and temperature at any given depth is higher for the previous model, when all other parameters are held constant. This also yield a smaller required depth for a desired temperature, although the difference becomes negligible as the surface heat flux becomes larger. The maximum differences of two models lay around \(q_s=50 \dfrac{mW}{m^2}\).
The solution of the PDE \[\dfrac{\partial T}{\partial t}= \dfrac{\lambda}{\rho c}\nabla ^{2}T=\alpha\nabla ^{2}T\] can be obtained by the Fourier Transform: \[\dfrac{\partial}{\partial t}[F_{z\rightarrow s}[T]]=-4\pi^2\alpha s^2 F[T]\] \[F[T](s,t)=F[T](s,0)e^{-4\pi^2\alpha s^2 t}\] \[T(z,t)=T(z,0)*F^{-1}[e^{-4\pi^2\alpha s^2 t}]\]
Where \(F^{-1}[e^{-4\pi\alpha s^2 t}]\) is \[F^{-1}[e^{-4\pi\alpha s^2 t}]=\int ^{\infty}_{-\infty} e^{-4\pi^2\alpha s^2 t}e^{i2\pi sz}ds\\=\int ^{\infty}_{-\infty} e^{-4\pi^2\alpha t(s+\dfrac{z}{i4\pi\alpha t})^2}e^{-\dfrac{z^2}{4\alpha t}}ds=\dfrac{1}{\sqrt{4\pi\alpha t}}e^{-\dfrac{z^2}{4\alpha t}}\]
Thus \[T(z,t)=T(z,0)*F^{-1}[e^{-4\pi^2\alpha s^2 t}]\\=\dfrac{1}{\sqrt{4\pi\alpha t}}\int ^{\infty}_{-\infty} T(z',0)e^{-\dfrac{(z-z')^2}{4\alpha t}}dz'\]
For our initial and boundary conditions we shall have \[T(z,t)\\=\dfrac{T_a}{\sqrt{4\pi\alpha t}}(\int ^{\infty}_{0}e^{-\dfrac{(z-z')^2}{4\alpha t}}dz'-\int ^{0}_{-\infty}e^{-\dfrac{(z-z')^2}{4\alpha t}}dz')\\=\dfrac{T_a}{\sqrt{4\pi\alpha t}}(-1+2\int ^{z}_{-\infty}e^{-\dfrac{(z-z')^2}{4\alpha t}}dz')\\=T_a(-1+2I(\dfrac{z}{\sqrt{2\alpha t}}))\]
Since R is a statistic software, it is easier to view the integral as a cumulative probability function of a normal distributed random variable with mean of zero and variance of \(2\alpha t\).
The resulting evolution of temperature profile is given below:
\[q_s(t)=-q(0,t)=\lambda\dfrac{\partial T}{\partial z}|_{z=0}=\dfrac{\lambda T_a}{\sqrt{\pi\alpha t}}\]
The characteristic length of the system can be obtained given the surface heat flux: \[L^*=\sqrt{\alpha t}=\dfrac{\lambda T_a}{\sqrt{\pi}q_s}\] From the graph on the next page, it is obvious that given same surface flux, oceanic geotherms have higher temperature profiles than those of continental geotherms at any depth z>1 km. This is because heat from Earth’s mantle is convected directly from the mid ocean ridge, a more efficient way of transfering heat and increasing temperature than by radiogenic sources.
To make the best use of the statistic functions of R, we rescale \(T_z=180\) to \(F=0.5+180/2600=0.5692\), the resulting quantile for a standardized normal distribution would be \(q=qnorm(F,0,1)=0.1744\).
As shown in 3(a), simply multiplying this by \(\sqrt{2\alpha t}\) yields the corresponding depth. The characteristic length \(\sqrt{\alpha t}\) is derived from 3(c).
As seen on the graph next page, the required depth to obtain a certain temperature is less for oceanic geotherms. The difference is also maximized around \(q_s=50 \dfrac{mW}{m^2}\)
| Month | Average.Temperature | HDD | Average.Power.per.Area |
|---|---|---|---|
| July | 19.3 | 4 | 0.0969060 |
| August | 18.6 | 27 | 0.6541154 |
| September | 14.8 | 127 | 3.1793239 |
| October | 11.0 | 266 | 6.4442480 |
| November | 6.6 | 438 | 10.9649123 |
| December | 3.1 | 493 | 11.9436627 |
| January | 1.4 | 474 | 11.4833593 |
| February | 0.5 | 390 | 10.4606551 |
| March | 5.8 | 304 | 7.3648549 |
| April | 10.4 | 168 | 4.2057198 |
| May | 13.4 | 63 | 1.5262693 |
| June | 17.4 | 20 | 0.5006809 |
Discussion: The average power per living area is negatively correlated to the average temeperature of that month. This means that a colder month leads to more energy demand, which is very common in high latitude regions.
The general form of the fitting problem is to solve: \[\Sigma_{i} (\widehat{y}_i-y_i) \dfrac{\partial \widehat{y}_i}{\partial \theta_j}=0\] For every fitting parameter \(\theta_j\).
In our particular case \(T(t)=a_1+a_2cos(\omega_y t)+a_3sin(\omega_y t)\), setting \(a_1=0\) (without loss of generality) and we have: \[\begin{cases} (\Sigma_{i}cos^2(\omega_y t_i))a_2+(\Sigma_{i}cos(\omega_y t_i)sin(\omega_y t_i))a_3=\Sigma_{i}cos(\omega_y t_i)y_i\\ (\Sigma_{i}cos(\omega_y t_i)sin(\omega_y t_i))a_2+(\Sigma_{i}sin^2(\omega_y t_i))a_3=\Sigma_{i}sin(\omega_y t_i)y_i\end{cases}\]
Since all the coefficients are constant, this is a simple linear system to be solved. Mind that \(t_i\) should be taken as the midpoint between two months (i.e. 0.5 month, 1.5 month, so on) to get the best fit.
Due to our constant sampling rate, the original linear system can be reduced into two independent equations: \[a_2=\dfrac{\Sigma_{i}cos(\omega_y t_i)y_i}{(\Sigma_{i}cos^2(\omega_y t_i))}\\a_3=\dfrac{\Sigma_{i}sin(\omega_y t_i)y_i}{(\Sigma_{i}sin^2(\omega_y t_i))}\]
Which yields \(a_2=8.4806\) and \(a_3=3.0325\). The alternative form \(T(t)=T_m+A_ycos(\omega(t-t_y))\) yields \(T_m=10.1917\), \(A_y=9.0065\) and \(t_y=0.6559\).
The general solution to the time dependent heat equation can be decomposed to a transcient part and a stationary part. In assignment 3 we have derived the solution to the transcient part. We now consider the stationary part of the general solution.
Suppose the stationary solution can be written as \[T_{st}(z,t)=\int^{\infty}_{-\infty}A(z,f)e^{i2\pi ft}df\]
Then the heat equation can be rearrange into \[\int^{\infty}_{-\infty}(i2\pi f-\alpha\dfrac{\partial^2 A}{\partial z^2})e^{i2\pi ft}df=0\]
A non-trivial solution exists when \[i2\pi f-\alpha\dfrac{\partial^2 A}{\partial z^2}=0\\A=c_+e^{(1+i)\sqrt{\dfrac{\pi f}{\alpha}}z}+c_-e^{-(1+i)\sqrt{\dfrac{\pi f}{\alpha}}z}\] But only the second term is physically feasible (for we must have \(A \rightarrow 0 \ as \ z \rightarrow \infty\)). Thus we have \[A(z,f)=A(0,f)e^{-\sqrt{\dfrac{\pi f}{\alpha}}z}e^{-i\sqrt{\dfrac{\pi f}{\alpha}}z}\\T_{st}(z,t)=\int^{\infty}_{-\infty}A(0,f)e^{-\dfrac{z}{d}}e^{i(2\pi ft-\dfrac{z}{d})}df\\d=\sqrt{\dfrac{\alpha}{\pi f}}\]
For our problem we have \[A(0,f)=A_ye^{-i2\pi ft_y}\delta(f-f_y)\]
So that \[T_{st}(z,t)=A_ye^{-\dfrac{z}{d}}e^{i(2\pi f_y(t-t_y)-\dfrac{z}{d})}\] The resulting subsurface temperature profiles are given below.