-
Notifications
You must be signed in to change notification settings - Fork 3
/
geotherms.m
190 lines (158 loc) · 5.66 KB
/
geotherms.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
function varargout=geotherms(Ts,qs,D,k,H)
% [Tz,z,qz]=GEOTHERMS(Ts,qs,D,k,H)
%
% Illustrates the behavior of some geotherms
%
% INPUT:
%
% Ts surface temperature [K]
% qs surface heatflow [W/m/m]
% D depth extent of the top medium [km]
% k thermal conductivity [W/m/K]
% H volumetric heat production [W/m/m/m]
%
% OUTPUT
%
% Tz temperature as a function of depth
% z the depth
% qz heat flow as a function of depth
% p plot handles
%
% Last modified by fjsimons-at-alum.mit.edu, 11/27/2021
% Surface temperature
defval('Ts',15);
% Surface heatflow
defval('qs',40e-3);
% Volumetric heat production
defval('H',5.9e-7);
% Thermal conductivity
defval('k',2.5);
% The depth extent of the top layer
defval('D',10000)
% Output units
uni='km';
if strcmp(uni,'km')
% The divisor that goes from input to output
unf=1000; else unf=1;
end
% Set up the query depths
Dmul=3;
defval('z1',linspace(0, D))
defval('z2',linspace(D,Dmul*D))
z=[z1(:); z2(:)];
% Now make the axis handles
ah=krijetem(subnum(2,3));
% An absolute reference for where the heat production will be plotted
% if left [] it will be relative to the maximum temperature reached
mTz=550;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FIRST CASE
% Calculate the first temperature profile
Tz1=Ts +qs /k*z1 -H/2/k*z1.^2;
% Calculate the first heat flow profile
qz1= qs -H*z1;
% Calculate the second temperature profile as a simple continuation
Tz2=Tz1(end)+qz1(end)/k*(z2-D)-H/2/k*(z2-D).^2;
% Calculate the second heat flow profile as a simple continuation
qz2= qz1(end) -H*(z2-D);
[Tz,qz,Tr,Trange,mTz]=assembly(Tz1,Tz2,qz1,qz2,mTz);
axes(ah(1))
p(1:4)=plotem(z,Tz,qz,unf,D,Trange,k,qz1);
p(5)=plot([mTz+Tr mTz+Tr NaN mTz+Tr mTz+Tr],...
[0 D NaN D z2(end)]/unf,'LineWidth',1.5,'Color','b');
hold off
[xl(1),yl(1),tl(1)]=...
labeling(ah(1),z,Tz,qz,Ts,qs,H,k,z1,Tz1,qz1,uni,unf,Trange);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SECOND CASE
% Calculate the first temperature profile
Tz1=Ts +qs /k*z1 -H/2/k*z1.^2;
% Calculate the first heat flow profile
qz1= qs -H*z1;
% Calculate the second temperature profile without any heat producers
Tz2=Tz1(end)+qz1(end)/k*(z2-D)-0/2/k*(z2-D).^2;
% Calculate the second heat flow profile - keep the zero for dimensions
qz2= qz1(end) -0*(z2-D);
[Tz,qz,Tr,Trange,mTz]=assembly(Tz1,Tz2,qz1,qz2,mTz);
axes(ah(2))
p(6:9)=plotem(z,Tz,qz,unf,D,Trange,k,qz1);
p(10)=plot([mTz+Tr mTz+Tr 0 0],...
[0 D D z2(end)]/unf,'LineWidth',1.5,'Color','b');
hold off
[xl(2),yl(2),tl(2)]=...
labeling(ah(2),z,Tz,qz,Ts,qs,H,k,z1,Tz1,qz1,uni,unf,Trange);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% THIRD CASE
% Calculate the first temperature profile
Tz1=Ts +(qs-D*H) /k*z1 +D^2*H/k.*(1-exp(-z1/D));
% Calculate the first heat flow profile
qz1= qs -D*H*(1-exp(-z1/D));
% Calculate the second temperature profile as a simple continuation
Tz2=Tz1(end)+(qz1(end)-D*H)/k*(z2-D)+D^2*H/k.*(1-exp(-(z2-D)/D));
% Calculate the second heat flow profile - keep the zero for dimensions
qz2= qz1(end) -D*H*(1-exp(-z1/D));
[Tz,qz,Tr,Trange,mTz]=assembly(Tz1,Tz2,qz1,qz2,mTz);
axes(ah(3))
p(11:14)=plotem(z,Tz,qz,unf,D,Trange,k,qz1);
p(15)=plot([mTz+Tr]*exp(-z/D),z/unf,'LineWidth',1.5,'Color','b');
hold off
[xl(3),yl(3),tl(3)]=...
labeling(ah(3),z,Tz,qz,Ts,qs,H,k,z1,Tz1,qz1,uni,unf,Trange);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Common Cosmetics
set(ah,'ytick',[0:Dmul]*D/unf)
set(ah,'xgrid','on')
longticks(ah,2)
delete(ah(4:end)); ah=ah(1:3);
movev(ah(1:3),-0.2)
for index=1:3
xel(index,:)=get(ah(index),'xlim');
end
set(ah,'xlim',[min(xel(:,1)) max(xel(:,2))]);
% set(ah(index),'position',get(ah(index),'position')+[0 -0.2 0 0])
% Optional output
varns={Tz,z,qz};
varargout=varns(1:nargout);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [Tz,qz,Tr,Trange,mTz]=assembly(Tz1,Tz2,qz1,qz2,mTz)
% Put them together
Tz=[Tz1(:) ; Tz2(:)];
qz=[qz1(:) ; qz2(:)];
% Reasonable offset
Tr=range(Tz)/10;
% Relative unless absolute input
defval('mTz',max(Tz));
% Scaling to render qz beautiful
if length(unique(qz(:)))~=1
% qz=scale([qz ; 0],[min(Tz)+Tr mTz-Tr]);
qz=scale([qz ; 0],[0 mTz-Tr]);
qz=qz(1:end-1);
else
qz=repmat(median(Tz),size(qz));
end
% Temperature axis
Trange=[-mTz/20 mTz+2*Tr];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function p=plotem(z,Tz,qz,unf,D,Trange,k,qz1)
p(1)=plot(Tz,z/unf,'LineWidth',1.5,'Color','r');
hold on
p(2)=plot(qz,z/unf,'LineWidth',1.5,'Color','g');
p(3)=plot(Trange,[D D]/unf,'--','LineWidth',0.5,'Color',grey);
% Plot a tangent - not that qz was already scaled, going in
%plot([Tz(1) Tz(1)+[Tz(2)-Tz(1)]/[z(2)-z(1)]*z(end)],[0 z(end)/unf])
p(4)=plot([Tz(1) Tz(1)+qz1(1)/k*z(end)],[0 z(end)/unf],'Color',grey);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [xl,yl,tl]=labeling(ah,z,Tz,qz,Ts,qs,H,k,z1,Tz1,qz1,uni,unf,Trange)
axis ij
ylim([0 max(z/unf)])
xlim(Trange)
set(ah,'xtick',sort(round([Ts Tz1(end) round(qz(length(z1))) max(Tz)])))
xl=ylabel(sprintf('depth (%s)',uni));
yl=xlabel('temperature | heat flow | heat production');
tl=title(sprintf(...
'%s %i%s/km\n%s %i mW/m^3\n%s %i mW/m^3\n%s %3.1f%sW/m^3\n',...
'surface thermal gradient',round(qs/k*unf),176,...
'surface heat flow',round(qs*1000),...
'reduced heat flow',round(qz1(end)*1000),...
'heat production',H*1e6,'\mu'),...
'FontWeight','normal');