% MF_ADAPT -- Pipe Cooling using Seawat for Viscosity and heat transport
%  TO 100405

clear variables
close all

basename='Geotherm';

% NOTICE in this model everything is in years !!!

% This exampple shows dissipation of the temperature anomaly in a geothermal
% heat system, caused by cooling the geothermal aquifer due to back injection
% of the cooled water. We start with an initial temperature distribution
% above and below the aquifer at the point considered. This point is
% characterized by the time since the cold front passed this point.
% We will compute the temperature dissipation using a single column of
% cells, starting with an initial temperature.

%% Parameters
L     =2000;    % [m] half depth of system
W     =100;     % [m] thickness of aquifer
Z0    =   0;    % [m] elevation of center of geothermal aquifer
T0    =  10;    % [C] temperature at ground surface
peff  = 0.20;   % [ ]  effective porosity
rhos  =2650;    % [kg/m3] density of solids
rhow  =1000;    % [kg/m3] density of water
rhob  =(1-peff)*rhos; % [kg/m3] bulk density of porous medium
k     = 40;     % [m/y] hydraulic conductivity
Lw    =0.6*30e6; % [J/y/m/K] heat conductance of water
Ls    =3.0*30e6; % [J/y/m/K] heat conductance of solides
cw    =4200;     % [J/kg/K]  heat capacity of the water
cs    = 800;     % [J/kg/K]  heat capacity of the solids
Lambda=peff*Lw+(1-peff)*Ls ;            % [J/y/m/K] bulk heat conductance
rhoc  =peff*rhow*cw +(1-peff)*rhob*cs;  % [J/m3/K]  bulk heat capacity
DH    =Lambda/(peff*rhow*cw);     % [m2/y]  diffusivity for heat
D     =DH;
Kd=cs/(rhow*cw);      % [m3/kg] distribution coefficient for heat
R =1+rhob*Kd/peff;    % [-] retardation
DTEMP=-50;            % [C] initial temperature
C0   =DTEMP;          % [C] initial temperature or concentration
C0max=DTEMP;          % [C] maximum temperature in plot of results
G    =30/1000;        % [K/m] geothermal temperature gradient
tp   =20;             % [y] time since front passage at end of project life

%% Mesh
yGr=[sinespace(-L,0,2000,pi/2,pi/20), sinespace(0,L,2000,pi/20,pi/2)];
zGr=[-1 1]/2;
xGr=[-1 1]/2;

[xGr,yGr,zGr,xm,ym,zm,Dx,Dy,Dz,NX,NY,NZ]=modelsize3(xGr,yGr,zGr);

Z=zGr(ones(NY,1),ones(1,NX),1:NZ+1);

[XM,YM,ZM]=meshgrid(xm,ym,ym);
[DX,DY,DZ]=meshgrid(Dx,Dy,Dz);

%% Generate all other matrices
IBOUND=ones( NY,NX,NZ);
ICBUND=ones( NY,NX,NZ);  ICBUND([1 end],1,1)=-1; % Ground Surface and z=-4000m
STRTHD=zeros(NY,NX,NZ)*0;
STCONC=zeros(NY,NX,NZ)*0;
HK    =ones(NY,NX,NZ)*k;
VK    =HK;
PEFF  =ones(NY,NX,NZ) *peff;
DMCOEF=ones(size(IBOUND))*DH;  % m2/y

%% Initial concentration
zTop=Z0+W/2;
zBot=Z0-W/2;

T=T0+G*(ym(1)-ym);

DT=zeros(size(ym));

sigma=sqrt(2*DH/R*tp);

DT(ym>=zTop)=(DTEMP+G*W/2)*erfc( (ym(ym>=zTop)-zTop)/(sigma*sqrt(2)));
DT(ym<=zBot)=(DTEMP-G*W/2)*erfc( (zBot-ym(ym<=zBot))/(sigma*sqrt(2)));
I=find(ym>=zBot & ym<=zTop);
DT(I)=interp1([ym(I(1)-1) ym(I(end)+1)],[DT(I(1)-1), DT(I(end)+1)],ym(I));

STCONC(:,:,:)=T+DT;

%plot(squeeze(ym),squeeze(STCONC))

%% Save variable that we need in mf_analyze.m
save underneath.mat G DTEMP T0 D R W L tp 

