function[] = cB_evol(scen,runnum,tstart,tend,dt,dts,GTfac,Dfac,zpmin,zpmax) % solution of cB evolution Eqn. (30) (supplemented by (4), (5), (9) and (14)) % on the spatial domain [zpmin to zpmax] at time step dt % with background fields of model ice-core parameters (i.e. age-depth profile and depth profiles of w, T and dg) % in two different core scenarios (scen): GRIP, EPIC % % runnum = Run Number % % GTfac (Gibbs-Thomson Factor): 0 turn off Gibbs-Thomson effect % 1 turn on Gibbs-Thomson effect % Dfac (Diffusivity Suppression Factor): 0 to 1 % (1) set up % parameters yr = 365.25*24*3600; D = 5e-10 * Dfac; % m2 s-1 vein impurity diffusivity rhoi = 917; rhow = 1e3; r = rhoi/rhow; % age axis set-up t = [tstart:dt:tend+dt]; % in years Lt = length(t); % time sampling of results to store in the output results file % dts is the time interval of sampling of the results ts = [tstart:dts:tend]; Lts = length(ts); % material reference frame (zp means z-prime) and associated model variables dzp = 0.0025; zp = [zpmin : dzp : zpmax]; N = length(zp); tapvec = tapervec(length(zp),100); % allocate empty matrices for key variables cB = zeros(1,N); cBn = zeros(1,N); % for current time step, and next time step cBout = zeros(Lts,N); cout = zeros(Lts,N); rvout = zeros(Lts,N); phiout = zeros(Lts,N); ageout = zeros(Lts,N); wcout = zeros(Lts,N); wrout = zeros(Lts,N); kapout = zeros(Lts,N); mout = zeros(Lts,N); qout = zeros(Lts,N); zout = zeros(1,Lts); % ice-core parameters and core reference frame eval(['load core_data_' scen]); % load z0-t0 scale and background fields for the model ice core cB0 = 1; % baseline bulk impurity concentration (micro M) in the simulation z = interp1(t0,z0,t); % depth of material at time/age t w = interp1(z0,w0,z); % ice velocity of material at depth z (at time/age t) % calculate baseline (steady state) vein network conditions [rv0,c0,phi0] = calc_vein_nonlin_1(dg0,T0,cB0,z0,GTfac); m0 = rhoi * w0 .* [phi0(2)-phi0(1) , 0.5*(phi0(3:end)-phi0(1:end-2)) , phi0(end)-phi0(end-1)]/dz; q0 = -((1-r)/rhoi) * dz * cumtrapz(m0); % (2) time evolution -- integrate the cB partial differential equation % initial condition for cB: apply the Gaussian doped peak (here, with amplitude 5 and width parameter 0.08) cB(1,:) = cB0 + 5*exp(-(zp/0.08).^2); for i = 1:Lt-1 % look up the local fields of temperature T, grain size dg, and relative velocity wr = w(g+zprime) - w(g) T = interp1(z0,T0,z(i)+[zp-0.5*dzp,zpmax+0.5*dzp]); % T on staggered grid dg = interp1(z0,dg0,z(i)+[zp-0.5*dzp,zpmax+0.5*dzp]); % dg on staggered grid wr = interp1(z0,w0,z(i)+zp) - w(i); % calculate vein-motion diffusivity, on staggered grid kap = calc_kappa(T); % calculate vein conditions, on staggered grid cBs = 0.5*[cB(N)+cB(1),cB(1:N-1)+cB(2:N),cB(N)+cB(1)]; [rvdum,cdum,phidum] = calc_vein_nonlin_1(dg,T,cBs,z(i)+[zp-0.5*dzp,zpmax+0.5*dzp],GTfac); % c = 0.5*(cdum(1:N)+cdum(2:N+1)); phi = 0.5*(phidum(1:N)+phidum(2:N+1)); % Rempel et al.'s velocity of 'anomalous diffusion', wc wc = -(1/dzp) * (D*yr) * (cdum(2:N+1)-cdum(1:N)) ./ c; % PDE % calculate advection terms on LHS by central difference cBprime1 = 0.5 * [ cB(2)-cB(N) , cB(3:N)-cB(1:N-2) , cB(1)-cB(N-1) ] / dzp; term1 = (r*wr + wc) .* cBprime1; % calculate residual diffusion (kappa) term on RHS cBprime2 = [ cB(1)-cB(N) , cB(2:N)-cB(1:N-1) , cB(1)-cB(N) ] / dzp; dum = kap(find(zp==0)) .* cBprime2; term2 = (dum(2:N+1)-dum(1:N)) / dzp; % calculate -cB*(dwc/dz) term on RHS wcprime = [wc(2)-wc(1), 0.5*(wc(3:N)-wc(1:N-2)),wc(N)-wc(N-1)] / dzp; term3 = - cB.*wcprime; % time stepping cBn = cB + (dt/r) * (term2 + term3 - term1) .* tapvec; % % sample and store results if rem((i-1),round(dts/dt)) == 0 j = 1 + (i-1)/round(dts/dt); zout(j) = z(i); cBout(j,:) = cB; cout(j,:) = c; rvout(j,:) = 0.5*(rvdum(1:N)+rvdum(2:N+1)); phiout(j,:) = phi; wcout(j,:) = wc; wrout(j,:) = wr; kapout(j,:) = 0.5*(kap(1:N)+kap(2:N+1)); ageout(j,:) = interp1(z,t,z(i)+zp); % % calculate m and q Tn = interp1(z0,T0,z(i+1)+zp); dgn = interp1(z0,dg0,z(i+1)+zp); [rvn,cn,phin] = calc_vein_nonlin_1(dgn,Tn,cBn,z(i+1)+zp,GTfac); dum = kap .* [ phi(1)-phi(N) , phi(2:N)-phi(1:N-1) , phi(1)-phi(N) ] / dzp; m = rhoi * ( (phin-phi)/dt + wr.*(phidum(2:N+1)-phidum(1:N))/dzp - (dum(2:N+1)-dum(1:N))/dzp ); mout(j,:) = m; qout(j,:) = interp1(z0,q0,z(i)+zpmin) - ((1-r)/rhoi)*dzp*cumtrapz(m); % % j % pause(0.001) end % move to next time step cB = cBn; % end clear i j T dg cB cBn cBs rvdum cdum phidum c phi wc wr kap term1 term2 term3 cBprime1 cBprime2 dum wcprime m Tn dgn rvn cn phin tapvec phi0 m0 q0 rv0 c0 cB0 clear yr D H h n a w Lt t z dz N rhoi rhow r Lts ageout kapout wrout eval(['save results_' scen '_' num2str(runnum)])