-------------------------------------------------------------------------------
"Pervasive diffusion of climate signals recorded in ice-vein ionic impurities"
by Felix S. L. Ng
-------------------------------------------------------------------------------

This ORDA respository contains the MATLAB computer code and the simulated data of 
key numerical experiments described in the above manuscript.

Reference: 
Ng, F. S. L. (2021)
Pervasive diffusion of climate signals recorded in ice-vein ionic impurities.
The Cryosphere Discussions, https://tc.copernicus.org/preprints/tc-2020-217/


License for the data files: CC BY 4.0.
License for the code: MIT open source license.


-------------------------------------------------------------------------------
(1) Outline

The MATLAB code supplied here solves Eqn. (30) of the paper, which describes the spatio-
temporal evolution of a signal in the bulk impurity concentration, cB, in an ice core. 
The spatial coordinate, zp (z' in the paper), is a material-following coordinate system
that follows a short section of ice as it descends the ice column towards the bed of the 
ice sheet. See the paper for more information.

The computer code consists of the following files:

Main program:

	cB_evol.m		integration of the partial differential equation for cB(z',t) (Eqn. (30))

Subsidiary programs:

	calc_kappa.m		calculate residual diffusivity, kappa

	calc_vein_nonlin_1.m	calculate vein conditions (impurity concentration c, vein-face radius of curvature rv) 
                              		and porosity phi, for ice at temperature T and depth z that
					carries the bulk impurity concentration cB and has grain size dg
				
	tapervec.m		subroutine used within cB_evol.m, to suppress problems at the domain boundaries (zpmin, zpmax)
					when a simulated signal has migrated near zpmin or zpmax, to interfere with 
					the zero-gradient boundary there

	test_runs.m		script file for recreating four key experiments described in the paper; see Item (2)

	plot_results_movie.m	script file for viewing movies of the simulated results for cB; see Item (2)


Input datafiles, containing the background glaciological fields at two model ice-core sites:

	core_data_GRIP.mat	GRIP core site background fields

	core_data_EPIC.mat	EPICA Dome C core site GRIP background fields

	(see further explanations in Section 3.1 and Fig. 4 of the paper).

	
Output datafiles, containing the computed results of four key experiments (MATLAB .mat format):

	results_GRIP_1.mat	GRIP control run (see Fig. 5a-c of paper)

	results_GRIP_2.mat	GRIP run with gamma = 0 (Fig. 5d-f)

	results_EPIC_3.mat	EPICA control run (Fig. 6a-c)

	results_EPIC_4.mat	EPICA run with gamma = 0 (Fig. 6d-f)


-------------------------------------------------------------------------------
(2) Running experiments and viewing their simulated results

To repeat the four key experiments, run the following command:

	>> test_runs;

To view movies of the computed time-evolving signal (cB) of the experiments in MATLAB, 
run the following commands:

	>> plot_results_movie('results_GRIP_1');

	>> plot_results_movie('results_GRIP_2');

	>> plot_results_movie('results_EPIC_3');

	>> plot_results_movie('results_EPIC_4');


-------------------------------------------------------------------------------
(3) Variables in the simulation and in the input and output datafiles

(a) Preliminaries:

	The spatial domain is defined by the vector
		zp = [zpmin: dzp: zpmax] ,
	where dzp is spatial step size. The default is dzp = 0.0025 m (line 32 of cB_evol.m).
	All distances are in metres.
	
	The time steps of the numerical integration are defined by the vector
		ts = [tstart: dt: tend] ,
	where dt is the time step. All times are in years.
	dt is specified when calling cB_evol.m. See test_runs.m, for example.
	
	The time steps sampled for the output variable fields are defined by
		ts = [tstart: dts: tend] ,
	where dts is the sampling time step (> dt).

	Each output variable field (e.g. cB) is presented as a matrix 
	with N rows and M columns, where N = length(ts) and M = length(zp).

(b) Input datafiles -- background fields at model ice-core site:

	Variables	Explanation / notes / units
	---------	---------------------------

	H		ice column thickness (m)

	a		surface accumulation rate (ice equivalent) (m yr^{-1}))

	h		kink height about the bed (m) [in the Dansgaard-Johnsen flow model; see Section 3.1 and Fig. 4 of the paper]

	n		exponent [in the flow model of Ritz (1992); see Section 3.1 and Fig. 4 of the paper ]

	m 		basal melt rate (m yr^{-1})) [in the flow model of Ritz (1992); see Section 3.1 and Fig. 4 of the paper]

	z0            	depth (m); this vector has the same size as t0, T0, dg0, w0

	dz		spatial step size in z0 (m)

	t0 		age of the ice (a, or yr in terms of the time of simulation) at depth z0 
	
	T0            	ice temperature (degC) at depth z0 

	dg0		mean grain size (mm) at depth z0 

	w0            	ice submergence velocity (m yr^{-1}) at depth z0 

(c) Output datafiles:

	Variables	Explanation / notes / units
	---------	---------------------------

	Simulation parameters:

	runnum        	run number. Decided by user when calling cB_evol.m; see test_runs.m, for example.

	scen          	= either 'GRIP' or 'EPIC'. This defines which input datafile to use 
			as the source of the ice-core background fields. Decided by user when 
			calling cB_evol.m; see test_runs.m, for example.

	Dfac          	suppression factor on the molecular diffusivity of impurity in vein water (D);
			Dfac is specified a value from 0 to 1. Dfac = 1 for no suppression.
			Decided by user when calling cB_evol.m; see test_runs.m, for example.

	GTfac         	factor controlling whether the Gibbs-Thomson effect is switched on 
			(GTfac = 1) or off (GTfac = 0). Decided by user when calling cB_evol.m; 
			see test_runs.m, for example.

	ts 		time vector (of the output samples); see item (3a) 

	tstart		initial time in the simulation (yr). This is the age of the ice where
			the doped signal is applied. Decided by user when calling cB_evol.m;
			see test_runs.m, for example.

	tend		end time in the simulation (yr). Decided by user when calling cB_evol.m;
			see test_runs.m, for example.

	dts		sampling time step (yr); see item (3a). Decided by user when calling cB_evol.m; 
			see test_runs.m, for example.

	dt 		timestep of integration in the simulation (yr). Decided by user when calling 
			cB_evol.m; see test_runs.m, for example.

	zp		space vector, defining the simulation domain (= z' in the paper); see item (3).

	zpmin		position of the left-hand edge of the simulation domain (m)
			Decided by user when calling cB_evol.m; see test_runs.m, for example.

	zpmax		position of the right-hand edge of the simulation domain (m)
			Decided by user when calling cB_evol.m; see test_runs.m, for example.

	dzp 		spatial step size (m). Decided by user when calling cB_evol.m;
			see test_runs.m, for example.

 	Vectors of background fields of the model ice-core site that have been used (e.g. for the purpose of
	interpolation) during the simulation:

	z0            	depth (m)

	t0 		age of the ice (a, or yr in terms of the time of simulation) at depth z0 
	
	T0            	ice temperature (degC) at depth z0 

	dg0		mean grain size (mm) at depth z0 

	w0            	ice submergence velocity (m yr^{-1}) at depth z0 

	Output variables: (see Item 3a)

  	zout        	depth of the ice section (m) (at zp = 0) at the time ts.
			(zout has the same length as ts.)

	cBout        	bulk impurity concentration (mu-M, i.e. 1e-6 mol per litre)

  	cout        	vein impurity concentration (M, i.e. mol per litre)

  	mout        	melt rate (kg m^{-3} yr^{-1})     
       
  	phiout      	porosity (dimensionless)

  	qout        	percolative water flux (m yr^{-1})

  	rvout       	vein-face radius of curvature (mu-m)

  	wcout       	Rempel et al. (2001) anomalous diffusion velocity; see Eqn. (14) of the paper.

-------------------------------------------------------------------------------