%% Values used in the paper 
% Natural reheaing of low-enthalpy geothermal aquifers
% TO 100816
clc

lambw=0.06;  fprintf('heat conductance of water =%g W/m/K\n',lambw);
lambs=3.0;   fprintf('heat conductance of solids=%g W/m/K\n',lambs);
eps = 0.2;   fprintf('porosity                  =%g [-]\n',eps);
lamb= eps*lambw+(1-eps)*lambs; % heat conductivity of bulkmatrix
             fprintf('bulk heat conductance     =%g W/m/K\n',lamb);
rhow= 1000;  fprintf('density of fluid (brine)  =%g kg/m3\n',rhow);
rhos= 2650;  fprintf('density of solids         =%g kg/m3\n',rhos);
cw  = 4200;  fprintf('heat capacity of water    =%g J/kg/K\n',cw);
cs  =  800;  fprintf('heat capacity of solids   =%g J/kg/K\n',cs);
rhowcw = rhow*cw;
             fprintf('vol heatcap of w rhowcw   =%g J/m3/K\n',rhowcw);   
rhoscs = rhos*cs;
             fprintf('vol heatcap of s rhoscs;  =%g J/m3/K\n',rhoscs);
rhoc   = eps*rhowcw+(1-eps)*rhos*cs;
             fprintf('vol heatcap of matrix rhoc=%g J/m3/K\n',rhoc);
G   = 0.03;  fprintf('geothermal gradient       =%g K/m\n',G);
Dtemp=50;    fprintf('temp drop in geoth. system=%g K\n',Dtemp);
W=100;       fprintf('thickness of aquifer      =%g m\n',W);
t = W*Dtemp*rhoc/(lamb*G); 
             fprintf(['time to heat aquifer from geothermal heat flow:\n',...
                 'in seconds %10g s\n',...
                 'in days    %10g d\n',...
                 'in years   %10g y\n'],t,t/365.25,t/(365.25*64600));
t_paper = W*Dtemp*round(rhoc/5e5)*5e5/(round(lamb*5)/5*G);
             fprintf(['time to heat aquifer from geothermal heat flow:\n',...
                 'in seconds %10g s\n',...
                 'in days    %10g d\n',...
                 'in years   %10g y\n'],t_paper,t_paper/365.25,t_paper/(365.25*64600));
             
DH = lamb/(eps*rhowcw); fprintf('DH = %g [m2/s] = % g [m2/y]\n',DH,DH*365.25*86400);
R  = rhoc/(eps*rhowcw); fprintf('Retardation    = % g [-]\n',R);
fprintf('DH/R = %g [m2/y]\n',DH*365.25*86400/R);
tt(1)=365.23*86400;  %   1 year
tt(2)=10*tt(1);      %  10 year
tt(3)=10*tt(2);      % 100 year
sigma=sqrt(2*DH/R*tt);
fprintf(['sigma (  1 y) = %g \n',...
         'sigma ( 10 y) = %g\n',...
         'sigma (100 y) = %g\n'],sigma);
t0=25*365.25*86400;
sigma0=sqrt(2*(DH/R)*t0); fprintf('sigma0 for t=%g y = %g [m]\n',t0/(365.25*86400),sigma0);
T0=50; T=0.05*T0;
sigma7=(T0/T)*(W+sigma0*sqrt(8/pi))/sqrt(2*pi);
                          fprintf('sigma(eq 7) = %g [m]\n,',sigma7);
t10=sigma7^2/(2*DH/R); fprintf('Reheating time= %g [y]\n',t10/(365.25*86400));
Q=365.25*24*200; fprintf('Q = %g m3/y\n',Q);
f=0.5;           fprintf('efficiency = %g [-]\n',f);
P=f*Q*Dtemp*rhowcw; fprintf('Power = %g J/yr = %g kWh/yr\n',P,P/(3600*1000));

             
