A. Code for the Algorithm of Fitzgibbon Et Al. (1999)

APPENDIX

A. Code for the algorithm of Fitzgibbon et al. (1999):

FitEllipse[data_]:=Module[{i=0,num=0,soln =0,D1={},D2={},S1={},S2={},S3={},C1={},T={},M={},a1={},a2={},a={},evals={},evects={}},

num=Length[data];

D1=Table[{data[[i]][[1]] data[[i]][[1]],data[[i]][[1]] data[[i]][[2]],data[[i]][[2]] data[[i]][[2]]},{i,1,num}];

D2=Table[{data[[i]][[1]],data[[i]][[2]],1},{i,1,num}];

S1=Transpose[D1].D1;

S2=Transpose[D1].D2;

S3=Transpose[D2].D2;

C1={{0,0,2},{0,-1,0},{2,0,0}};

T=-Inverse[S3]. Transpose[S2];

M=Inverse[C1].(S1+S2.T);

evals=Eigenvalues[M];

soln = Position[evals,Select[evals,#>0&][[1]]][[1,1]];

evects=Eigenvectors[M];

a1=evects[[soln]];

a2=T.a1;

a={a1[[1]],a1[[2]],a1[[3]],a2[[1]],a2[[2]],a2[[3]]};

a

]

B. Fluid mechanical derivations

The ‘Poisson equation’ of rectilinear flow of a Newtonian viscous fluid in the z-direction through an infinitely long parallel-wall inclined channel is given by (eqn 6.190 of Papanastasiou et al. 2000):

(∂2Uz/∂x2) + (∂2Uz/∂y2) =µ-1 [∂P/∂z - d1 g Sinθ] (1)

Here ‘x’ and ‘y’: are the mutually perpendicular directions that lie on the cross-section of the channel; ‘µ’: dynamic viscosity of the fluid; (∂P/∂x) is the pressure gradient along the channel leading to the extrusion of the fluid; ‘d1’: density of the fluid, ‘g’: acceleration due to gravity, and ‘θ’: inclination of the channel.

We now consider (i) the channel is very long but of finite length; (ii) the fluid rises up the channel due to the presence of another fluid (of density d0, > d1,) present beneath it; and (iii) the cross-section of the channel to be elliptical with ‘x’ and ‘y’ as the major- and the minor axes, and of lengths ‘2a’ and ‘2b’, respectively. The second consideration leads to ∂P/∂z = d g Sinθ

Therefore eqn (1) becomes

(∂2Uz/∂x2) + (∂2Uz/∂y2) = µ-1 g (d0 - d1) Sinθ (2)

Let Uz (x,y) be the velocity of the extruding fluid at the coordinate (x,y). Considering the wall of the channel to be static during fluid flow, the boundary condition is

Uz (x,y) = 0 at (x2 a-2 + y2 b-2) = 1 (3)

A dependent variable Uz/ is now introduced such that

Uz (x,y) = Uz/ (x,y) + c1 x2 + c2 y2 (4)

c1, c2 are constants, ≠ 0, and are to be solved so that the following conditions are satisfied: Uz/ (x,y) satisfies the Laplace equation, and (ii) Uz/ (x,y) is constant on the wall at a particular instant ‘t’.

Substituting eqn (4) into eqn (2),

(∂2Uz/∂x2) + (∂2Uz//∂y2) + 2 (c1 + c2) = µ-1 g (d0 - d1) Sinθ (5)

Uz/(x) will satisfy the Laplace Equation:

(∂2Uz//∂x2)+ (∂2 Uz//∂y2) = 0 (6)

if 2 (c1 + c2) = µ-1 g (d0 - d1) Sinθ (7)

From the boundary condition (eqn 3),

Uz/ (x,y)=-(c1 x2 + c2 y2) =-c1 (x2 + c2 c1-1 y2) ; at (x2 a-2 + y2 b-2) = 1 (8)

Fixing (c2 c1-1) = (a2 b-2) (9)

Uz/ (x,y) is constant at the channel boundary at any particular instant:

Uz/ (x,y)=-c1 a2 on (x2 a-2+y2 b-2) = 1 (10)

According to the maximum principle for the Laplace equation, Uz/ (x,y) has both its minimum and maximum values on the boundary of the domain. This means that Uz/ (x,y) is constant over the whole domain at a particular time:

Uz/ (x,y) = - c1 a2 (11)

Putting eqn (11) in eqn (4), and using eqn (9)

Uz (x,y) = (- c1 a2 + c1 x2 + c2 y2) = - c1 a2 (1 - x2 a-2 - c2 c1-1 y2 a-2) (12)

or, Uz (x,y) = - c1 a2 (1 - x2 a-2- y2 b-2) (13)

The constant c1 is obtained from eqn (7) and eqn (9)

c1 = 0.5 b2 µ-1 g Sinθ (d0 - d1) (a2 + b2)-1 (14)

Putting the c1 value of eqn (14) into eqn (13)

Uz (x,y) = -0.5 µ-1 g a2 b2 Sinθ (d0 - d1) (a2 + b2)-1 (1 - x2 a-2 - y2 b-2) (15)

Therefore, at the centre of extrusion (x = y= 0):

Uz (0,0) = 0.5 µ-1 g a2 b2 Sinθ (d0 - d1) (a2 + b2)-1 (16)

We now consider that an overburden of three immiscible layers of fluids with their lengths along the channel ‘h2’, ‘h3’ and ‘h4’ (see Fig. 6) and densities ‘d2’, ‘d3’ and ‘d4’ at the top of the fluid column of density ‘d1’. In such a case, eqn (1) is rewritten as:

(∂2Uz/∂x2) + (∂2Uz/∂y2) = µ-1 g Sinθ (d0 - d1 h1 H-1 – d2 h2 H-1 – d3 h3 H-1 – d4 h4 H-1) (17)

Here ‘H’ stands for the total thickness of the four layers (i.e. H = h1 + h2 + h3 + h4) and ‘h1’ is the thickness of the proto-TMC. Thus eqn (16) takes the following form:

Uz (0,0) = 0.5 µ-1 g a2 b2 Sinθ (d0 – d1 h1 H-1 - d2 h2 H-1 - d3 h3 H-1 - d4 h4 H-1) (a2 + b2)-1

(18)

The parameters ‘g’, ‘a’, ‘b’, ‘θ’, ‘d0’, ‘d1’, ‘d2’, ‘d3’, ‘d4’ and ‘H’ remaining constant,

Uz (0,0) ∞ µ-1 (19)

And, for the parameters ‘Uz (0,0)’ ‘g’, ‘a’, ‘b’, ‘d0’, ‘d1’, ‘d2’, ‘d3’, ‘d4’ and ‘H’ remaining constant,

µ ∞ θ (20)

C. Interval function of Mathematica used in calculating viscosity

CalcViscosity[listin_, uplift_] :=

Module[{i, h1, h2, h3, h4, H, theta, d, d1, d2, d3, d4, g = 980, mu,

uz = 0, a = 4500000, b = 1600000, md, md1, md2, md3, md4, ans,

anslist = {}, str0, strmin, strmax, mumin, mumax, reslist = {}},

(*Expecting data in the form h1,h2,h3,h4,theta,d,d1,d2,d3,d4 *)

uz = uplift/(10 365 24 60 60) ;

Print["Uplift rate = ", uplift, " mm/year, ", uz // N, " cm/s"];

For[i = 1, i <= Length[listin], i++,

(*Put data in correct format for analysis*)

h1 = listin[[i, 1]] 100000;

h2 = listin[[i, 2]] 100000;

h3 = listin[[i, 3]] 100000;

h4 = listin[[i, 4]] 100000;

theta = listin[[i, 5]];

d = listin[[i, 6]];

d1 = listin[[i, 7]];

d2 = listin[[i, 8]];

d3 = listin[[i, 9]];

d4 = listin[[i, 10]];

H = h1 + h2 + h3 + h4;

h1 = h1/Sin[theta Pi/180];

h2 = h2/Sin[theta Pi/180];

h3 = h3/Sin[theta Pi/180];

h4 = h4/Sin[theta Pi/180];

H = H/Sin[theta Pi/180];

(*Calculate viscosity*)

mu = 0.5 g a^2 b^2 Sin[

theta Pi/180] (d - d1 h1/H - d2 h2/H - d3 h3/H -

d4 h4/H)/((a^2 + b^2) uz);

(*Check Min and Max viscosity densities *)

mumax =

0.5 g a^2 b^2 Sin[

theta Pi/180] (Max[d] - Min[d1] h1/H - Min[d2] h2/H -

Min[d3] h3/H - Min[d4] h4/H)/((a^2 + b^2) uz);

mumin =

0.5 g a^2 b^2 Sin[

theta Pi/180] (Min[d] - Max[d1] h1/H - Max[d2] h2/H -

Max[d3] h3/H - Max[d4] h4/H)/((a^2 + b^2) uz);

(*output results*)

str0 =

"theta = " > ToString[theta] > ", h1 = " >

ToString[h1/100000 // N] > ", h2 = " >

ToString[h2/100000 // N] > ", h3 = " >

ToString[h3/100000 // N] > ", h4 = " >

ToString[h4/100000 // N] > ", H = " > ToString[H/100000 // N] ;

Print["***********************************************************************************\

"];

Print[str0];

Print["viscosity = ", mu];

Print["viscosity min = ", mumin];

Print["viscosity max = ", mumax];

(*Find a set of parameters for which the viscosity is zero*)

(*Calculate mean densities*)

md1 = Mean[d1[[1]]]; md2 = Mean[d2[[1]]]; md3 = Mean[d3[[1]]];

md4 = Mean[d4[[1]]];

(*find corresponding d which makes vicosity zero*)

ans = Solve[dx - md1 h1/H - md2 h2/H - md3 h3/H - md4 h4/H == 0,

dx][[1, 1, 2]];

(*output answer*)

Print[

"Example values of densities which make viscosity equal zero:"];

str0 =

"d = " > ToString[ans] > ", d1 = " > ToString[md1] >

", d2 = " > ToString[md2] > ", d3 = " > ToString[md3] >

", d4 = " > ToString[md4];

Print[str0];

str0 =

ToString[ans] > ", " > ToString[md1] > ", " > ToString[md2] >

", " > ToString[md3] > ", " > ToString[md4];

strmin =

ToString[Min[d]] > ", " > ToString[Max[d1]] > ", " >

ToString[Max[d2]] > ", " > ToString[Max[d3]] > ", " >

ToString[Max[d4]];

strmax =

ToString[Max[d]] > ", " > ToString[Min[d1]] > ", " >

ToString[Min[d2]] > ", " > ToString[Min[d3]] > ", " >

ToString[Min[d4]];

(* Assemble results for output for conversion to table form *)

AppendTo[

reslist, {theta, h1/100000, h2/100000, h3/100000, h4/100000,

H/100000, mumax, strmax, mumin, strmin, str0} // N];

];

reslist

]

1