About the Code: It Is Written in Quickbasic

ABOUT THE CODE: IT IS WRITTEN IN QUICKBASIC

It is important to note that the code does do – and what we did within the paper- is just to solve the constitutive equations with appropriate initial and boundary conditions. Thinking in terms of practical cases, the initial analysis of the bulk density profile (e.g. the fit to an exponential saturation function) is done with any standard software (Origin and even Excel can be used for this purpose). These fitting parameters can then be stated as initial conditions in the numerical code. The use of time-dependent sedimentation rates has to be introduced by appropriately rewriting a line in the code (several options are incorporated as comments in the code). The code does not solve the CRS model, but it can be adapted to read time series of sedimentation rates (HISTOW.DAT in the example) to be used them in calculations. The model output are two files (with a multicolumn structure) which can be plotted using conventional software: ZDK.DAT contains the final depth distribution of bulk density and the conductivity function, while TVW.DAT stores time series of sedimentation velocity and the sedimentation rate at SWI.

Useful links:

REM sentences are comments, as those after the symbol “ ‘ “

CLS

REM *********************************************************************

PRINT "THIS CODE SOLVES THE COINTINUITY EQUATION FOR BULK DENSITY (Eq.2)"

PRINT "UNDER SPATIAL GRADIENTS OF THE COMPACTION POTENTIAL GIVEN BY THE BUOYANCY"

PRINT "AND BEING THE SPATIAL GRADIENTS OF CONDUCTIVITY ROPORTIONAL TO"

PRINT "THE DEFAULT BULK DENSITY (Eqs. 5 AND 8)"

PRINT "Please, read the readme.doc file before running this code"

REM ---- THE CODE WORKS WITH ACTUAL DEPTH VARIABLE---

REM INITIAL CONDITIONS WILL BE DEFINED BY AN EXPONENTIAL-SATURATION FUNCTION

REM RHO(Z)= RHHO(INF)-RHO1*EXP(-ALPHA Z )

REM The parameter values have to be provided by the user

REM

REM ----BOUNDARY CONDITIONS AT THE SEDIMENT-WATER INTERFACE:

REM Sedimentation rate has to be provided with three options

REM a) A constant value b) Time series stored in a two column file

REMc) An analytical function of time (this has to be written within the code)

REM

REM TO COMPUTE CONDUCTITY THE SITE-SPECIFIC PROPORTIONALITY CONSTANT Cz

REM has to be provided by the user

REM

REMDEFAULT VALUES FOR SPATIAL AND TIME DISCRETIZATION ARE ASUMED IN THE CODE

REM but this can be changed by the user. Note that they have to fulfill

REM the first criterion for numerical stability.

REM THE TIME OF SIMULATION has to be provided by the user

REM *************************************************************************

REM *************************"opening output and input files"

OPEN "ZDK.dat" FOR OUTPUT AS #1

REM A 3-column file for output with z, rho(z), K(z) for t=t simulation

OPEN "TVW.dat" FOR OUTPUT AS #2

REM A 3-column file for output with t, v(t) -sed. velocity-, w(t)

INPUT "w(t) from file ? (1)=no, else =yes"; cond

IF cond > 1 THEN GOTO 10

GOTO 20

10 INPUT "name of the input file with t, w(t)? (eg. histow.dat)"; a$

REM expected format is "age [t], w(age)[g/cm^2y]"

OPEN a$ FOR INPUT AS #3

INPUT "length of data (e.g. 45)"; ld ' number of data points in the file

DIM w(ld, 2) 'Matrix storing ages and sedimentation rates

FOR j = ld TO 1 STEP -1

INPUT #3, a

w(j, 1) = a

INPUT #3, b

w(j, 2) = b

PRINT w(j, 1), w(j, 2)

NEXT j

CLOSE #3

GOTO 30

REM ***** dimensioning arrays and defining variables *************

REM Default dimension will be 300 (number of sections defined within the core)

20 PRINT "Please introduce constant w [g/cm^2y], for other options rewrite line 110"

INPUT "w? (e.g. 0.1)"; w0

30 DIM kz(300) ' This is the conductivity matrix

DIM d(300) ' Bulk density in actual time t

DIM dn(300) ' Bulk density in t+dt

dz = .1 ' Depth interval for discretization (can be changed by user)

PRINT "INTRODUCE BULK DENSITY PARAMETERS FOR INITIAL CONDITIONS"

INPUT "RHO(INF) [g/cm^3] (e.g. 0.46)"; rhoinf

INPUT "RHO-1 [g/cm^3] (e.g. 0.28)"; rho1

INPUT "SCALING FACTOR, alpha [1/cm] (e.g. 0.154)"; alfa

rho0 = rhoinf - rho1

dsol = 2.5 'Default value for density of solids (can be changed by user)

dw = 1 ' Default value for the density of water (can be changed by user)

beta = (dsol - dw) / dsol' Buoyancy

dt = .005 'Default value for time-step (can be changed by user)

INPUT "TIME OF SIMULATION [y] (e.g. 150)"; ts

lt = ts / dt

t = -dt

INPUT "VALUE OF Cz [1/y] FOR CONDUCTIVITY (e.g. 0.056)"; C

REM "INITAL CONDITIONS FOR BULK DENSITY (evaluated at mid-point interval)"

FOR s = 1 TO 300

Z = Z + dz

d(s) = rhoinf - rho1 * EXP(-alfa * (Z - dz / 2))

kz(300) = C / alfa * rho1 * EXP(-alfa * 300 * dz)' Value of K at largest depth

REM "computing initial values for conductivity from down to top"

FOR s = 299 TO 0 STEP -1

kz(s) = kz(s + 1) + C * dz * (rhoinf - d(s + 1))

REM ***** START THE MAIN BUCLE IN TME *****

IF cond > 1 THEN LET ij = 2: lij = w(ij, 1) ' This is for reading w from a file

FOR it = 1 TO lt

t = t + dt ' time is updated

IF cond = 1 THEN GOTO 110 'when w is not provided by a file

IF t < lij THEN GOTO 100

ij = ij + 1

lij = w(ij, 1)

100 w = w(ij - 1, 2) + (w(ij, 2) - w(ij - 1, 2)) / (w(ij, 1) - w(ij - 1, 1)) * (t - w(ij - 1, 1))

GOTO 150 ' To skip out the lines with other alternative definitions of w

REM ****** Alternative definitions for w *************

REM The one activated is the constant values"

REM To activate others option, comment with "REM" line 110 and write your function

REM Some examples are listed below after REM comments

110 w = w0 ' this is the constant value introduced by user

REM w = .1 + .05 * COS(2 * 3.14156 * t / 7) ' A periodic variation in w

REM w = .1 + .08 / (SQR(2 * 3.14159) * .2) * EXP(-(t - 140) ^ 2 / (2 * .2 ^ 2))

REM The previous was a Gaussian pulse

REM IF t > 130 THEN LET w = w + dt * .0006 'Activated with line 110 produces acceleration

REM ******* UPPER BOUNDARY CONDITION******************

REM We solve the sedimentation velocity v [cm/y]

150 v = (w - beta * kz(1)) / ((d(1) + d(2)) / 2)

dn(1) = d(1) ' updating bulk density

IF ABS(it / 10 - INT(it / 10)) = 0 THEN PRINT #2, t, v, w 'producing output

FOR s = 2 TO 299 ' ***** main bucle in actual depth variable

IF v >= 0 THEN LET dn(s) = d(s) + dt * (-v * ((d(s) + d(s + 1)) / 2 - (d(s) + d(s - 1)) / 2) - beta * (kz(s) - kz(s - 1))) / dz

IF v < 0 THEN LET dn(s) = d(s) + dt * (v * ((d(s) + d(s + 1)) / 2 - (d(s) + d(s - 1)) / 2) - beta * (kz(s) - kz(s - 1))) / dz

REM addapting upstream method / Check for negative dn(s) values not activated

dn(300) = dn(299) * 1! 'Boundary condition at largest depth

REM "This boundary condition applies when the stated lengt of the core is large enough

FOR s = 1 TO 300

d(s) = dn(s) 'new values are stored as old values for the next time step

REM "The conductivity matrix is updated for the new time step"

kz(300) = C / alfa * rho1 * EXP(-alfa * 300 * dz)

FOR s = 299 TO 0 STEP -1

kz(s) = kz(s + 1) + C * dz * (rhoinf - d(s + 1))

NEXT it

REM ****this is the end of the main bucle for time *******

1000 FOR s = 1 TO 300

PRINT #1, s * dz, d(s), kz(s) 'This produces output

CLOSE #1: CLOSE #2

PRINT "ZDK.dat and TVW.dat output files have been created. Bye!"

END