Comment on “Channel flow, tectonic overpressure, and exhumation of high-pressure rocks in the Greater Himalayas” by Marques et al. (2018)

Platt, John P.

doi:https://doi.org/10.5194/se-10-357-2019

Articles | Volume 10, issue 1

https://doi.org/10.5194/se-10-357-2019

© Author(s) 2019. This work is distributed under
the Creative Commons Attribution 4.0 License.

https://doi.org/10.5194/se-10-357-2019

© Author(s) 2019. This work is distributed under
the Creative Commons Attribution 4.0 License.

Articles | Volume 10, issue 1

Peer-reviewed comment

|

27 Feb 2019

Peer-reviewed comment |

| 27 Feb 2019

Comment on “Channel flow, tectonic overpressure, and exhumation of high-pressure rocks in the Greater Himalayas” by Marques et al. (2018)

John P. Platt

Download

Final revised paper (published on 27 Feb 2019)
Preprint (discussion started on 26 Sep 2018)
Corresponding article

Interactive discussion

Status: closed

AC: Author comment | RC: Referee comment | SC: Short comment | EC: Editor comment

- Printer-friendly version

- Supplement

RC1: 'Review of Comment by J.P. Platt on Marques et al. (2018), Channel flow, tectonic overpressure, and exhumation of high-pressure rocks in the Greater Himalayas', Stefan Markus Schmalholz, 25 Oct 2018
- SC1: 'Reply to Platt's comment', Fernando Ornelas Marques, 27 Oct 2018
  - AC2: 'Response to Reply by F.O Marques et al.', John Platt, 05 Nov 2018
- AC1: 'Response to Review by S.M. Schmalholz', John Platt, 29 Oct 2018
  - RC3: 'Reply to Platt', Stefan Markus Schmalholz, 26 Nov 2018
    - AC4: 'Response to Comment by Stefan Schmalholz', John Platt, 28 Nov 2018
RC2: 'Review of “Comment on Marques et al. (2018), Channel flow, tectonic overpressure, and exhumation of high-pressure rocks in the Greater Himalayas” by John P. Platt.', Evangelos Moulas, 21 Nov 2018
- AC3: 'Response to Review by E. Moulas', John Platt, 28 Nov 2018
  - RC4: 'Comment to “Response to Review by E. Moulas by J. P. Platt” by E. Moulas', Evangelos Moulas, 29 Nov 2018
    - AC5: 'Reply to Mousas', John Platt, 30 Nov 2018

Peer-review completion

AR: Author's response | RR: Referee report | ED: Editor decision

AR by John Platt on behalf of the Authors (04 Dec 2018) Author's response

ED: Reconsider after major revisions (11 Dec 2018) by Taras Gerya

Dear John,
In order to check you arguments about boundary conditions numerically (as suggested by one of the reviewers), I did numerical experiments for a weaker channel surrounded by stronger rocks under condition of free lithospheric surface present in the model with inclined gravity orthogonal to this surface. Rock rheology is viscoelastic. The respective matlab Viscoelastic2D_m code is copied below for your reference. GPa-level overpressure indeed forms in the channel at high channel viscosity (10^22 Pa*s) that are balanced by GPa-level deviatoric stresses (and under/over pressures) in surrounding stronger rocks. It remains questionable if such stresses can be sustained by these rocks but within the referred visco-elastic approximation the results are rather consistent with the paper of Marques.
Therefore I think that your statement that "The fixed upper boundary condition in their models violates force balance and is unphysical" is not correct. Having this boundary non-fixed does not preclude overpressure (if surrounding rocks are much stronger than the channel).
Therefore, I think your comment cannot be published in its present state.
With my best wishes.
Taras

% Solving Stokes and continuity eqs.
% in primitive variable formulation
% with variable viscosity
% for visco-elastic material
% using FD with staggered grid

% Clearing memory and figures
clear all; clf

% Define Numerical model
xsize=1000000; % Horizontal model size, m
ysize=300000; % Vertical model size, m
Nx=101; % Horizontal grid resolution
Ny=101; % Vertical grid resolution
Nx1=Nx+1;
Ny1=Ny+1;
dx=xsize/(Nx-1); % Horizontal grid step, m
dy=ysize/(Ny-1); % Vertical grid step, m

% Coordinates of different nodal points
% Basic nodes
x=0:dx:xsize; % Horizontal coordinates of basic grid points, m
y=0:dy:ysize; % Vertical coordinates of basic grid points, m
% Vx-Nodes
xvx=0:dx:xsize+dy; % Horizontal coordinates of vx grid points, m
yvx=-dy/2:dy:ysize+dy/2; % Vertical coordinates of vx grid points, m
% Vy-nodes
xvy=-dx/2:dx:xsize+dx/2; % Horizontal coordinates of vy grid points, m
yvy=0:dy:ysize+dy; % Vertical coordinates of vy grid points, m
% P-Nodes
xp=-dx/2:dx:xsize+dx/2; % Horizontal coordinates of P grid points, m
yp=-dy/2:dy:ysize+dy/2; % Vertical coordinates of P grid points, m

% Nodal arrays
% Basic nodes
ETA=zeros(Ny,Nx); % Viscosity, Pa*s
GGG=zeros(Ny,Nx); % Shear modulus, Pa
EXY=zeros(Ny,Nx); % EPSILONxy, 1/s
SXY=zeros(Ny,Nx); % SIGMAxy, 1/s
SXY0=zeros(Ny,Nx); % SIGMA0xy, 1/s
wyx=zeros(Ny,Nx); % Rotation rate, 1/s
% Vx-Nodes
RHOX=zeros(Ny1,Nx1); % Density, kg/m^3
KX=zeros(Ny1,Nx1); % Thermal conductivity, W/m/K
vx=zeros(Ny1,Nx1); % vx-velocity m/s
% Vy-Nodes
RHOY=zeros(Ny1,Nx1); % Density, kg/m^3
KY=zeros(Ny1,Nx1); % Thermal conductivity, W/m/K
vy=zeros(Ny1,Nx1); % vy-velocity m/s
% P-nodes
RHO=zeros(Ny1,Nx1); % Density, kg/m^3
RHOCP=zeros(Ny1,Nx1); % Volumetric heat capacity, J/m^3/K
ALPHA=zeros(Ny1,Nx1); % Thermal expansion, J/m^3/K
HR=zeros(Ny1,Nx1); % Radioactive heating, W/m^3
HA=zeros(Ny1,Nx1); % Adiabatic heating, W/m^3
HS=zeros(Ny1,Nx1); % Shear heating, W/m^3
ETAP=zeros(Ny1,Nx1); % Viscosity, Pa*s
GGGP=zeros(Ny1,Nx1); % Shear modulus, Pa
EXX=zeros(Ny,Nx); % EPSILONxx, 1/s
SXX=zeros(Ny,Nx); % SIGMA'xx, 1/s
SXX0=zeros(Ny,Nx); % SIGMA0'xx, 1/s
tk1=zeros(Ny1,Nx1); % Old temperature, K
tk2=zeros(Ny1,Nx1); % New temperature, K
vxp=zeros(Ny1,Nx1); % Vx in pressure nodes, m/s
vyp=zeros(Ny1,Nx1); % Vy in pressure nodes, m/s
pr=zeros(Ny1,Nx1); % Pressure, Pa

% Define markers
Nxm=(Nx-1)*4; % Marker grid resolution in horizontal direction
Nym=(Ny-1)*4; % Marker grid resolution in vertical direction
dxm=xsize/Nxm; % Marker grid step in horizontal direction,m
dym=ysize/Nym; % Marker grid step in vertical direction,m
marknum=Nxm*Nym; % Number of markers
xm=zeros(1,marknum); % Horizontal coordinates, m
ym=zeros(1,marknum); % Vertical coordinates, m
tm=zeros(1,marknum); % Material type
tkm=zeros(1,marknum); % Marker temperature, K
sxxm=zeros(1,marknum); % SIGMA'xx, Pa
sxym=zeros(1,marknum); % SIGMAxy, Pa

% Define properties of materials:
% lith. channel air
rhom = [3000 3000 1 ]; % Density, kg/m^3
etam = [1e+26 1e+22 1e+17 ]; % Viscosity, Pa s
gggm = [5e+10 5e+10 5e+10 ]; % Shear Modulus, Pa

% Model geometry
aplate=15; % plate dipping angle, deg
abase=30; % channel base angle, deg
atopo=0; % topography angle, deg
mchan=100000; % channel mouth width, m
bchan=150000; % channel base width, m
lchan=600000; % channel length (at the footwall), m
gx=9.81*sin(aplate/180*pi); % Horizontal gravity acceleration, m/s^2
gy=9.81*cos(aplate/180*pi); % Vertical gravity acceleration, m/s^2
% Channel hanging wall
xleft=mchan*cos((aplate+atopo)/180*pi);
yleft=ysize-mchan*sin((aplate+atopo)/180*pi);
xright=lchan+bchan*cos((abase)/180*pi);
yright=ysize-bchan*sin((abase)/180*pi);
% Lithostatic pressure
prl=zeros(Ny+1,Nx+1);
for j=1:1:Nx+1
for i=1:1:Ny+1
% Air/Lithosphere/channel
if(yp(i) % Air
prl(i,j)=1e+5;
else
% Surface interception point along gravity vector
y1=((yp(i)*tan(aplate/180*pi)-xp(j))*tan((aplate+atopo)/180*pi)+ysize)/...
(1+tan(aplate/180*pi)*tan((aplate+atopo)/180*pi));
x1=(ysize-y1)/tan((aplate+atopo)/180*pi);
% Distance to the surface along gravity vector
d=((xp(j)-x1)^2+(yp(i)-y1)^2)^0.5;
% Lithostatic pressire
prl(i,j)=1e+5+d*rhom(1)*9.81;
end
end
end

% Define marker coordinates, temperature and material type
rp=100000; % Plume radius, m
m=1; % Marker counter
for jm=1:1:Nxm
for im=1:1:Nym
% Define marker coordinates
xm(m)=dxm/2+(jm-1)*dxm+(rand-0.5)*dxm;
ym(m)=dym/2+(im-1)*dym+(rand-0.5)*dym;
% Marker properties
% Air/Lithosphere/channel
if(ym(m)<=ysize-xm(m)*tan((aplate+atopo)/180*pi))
% Air
tm(m)=3;
elseif(ym(m)>yleft+(xm(m)-xleft)/(xright-xleft)*(yright-yleft) &&...
xm(m) % Channel
tm(m)=2;
else
% Lithosphere
tm(m)=1;
end
% Update marker counter
m=m+1;
end
end

% Introducing scaled pressure
pscale=1e+21/dx;

% Define global matrixes
% Mechanical solution: L(), R()
N=Nx1*Ny1*3; % Global number of unknowns
L=sparse(N,N); % Matrix of coefficients (left part)
R=zeros(N,1); % Vector of right parts

% Mechanical boundary conditions: free slip=-1; No Slip=1
bcleft=-1;
bcright=-1;
bctop=-1;
bcbottom=1;
vplate=1.5e-9; % Plate velocity, m/s
jplate=fix(lchan/dx)-1; % max index for velocity BC

% Timestepping
dt=1e+12; % Computational timestep
dxymax=0.5; % Max marker movement per time step, grid steps
vpratio=1/3; % Weight of averaged velocity for moving markers
dsubgrids=1; % Subgrid stress diffusion parameter
timesum=0; % Time sum, s
for timestep=1:1:1000

% Interpolate properties from markers to nodes
% Basic nodes
ETASUM=zeros(Ny,Nx);
GGGSUM=zeros(Ny,Nx);
SXYSUM=zeros(Ny,Nx);
WTSUM=zeros(Ny,Nx);
% Vx-nodes
RHOXSUM=zeros(Ny1,Nx1);
WTXSUM=zeros(Ny1,Nx1);
% Vy-nodes
RHOYSUM=zeros(Ny1,Nx1);
WTYSUM=zeros(Ny1,Nx1);
% P-Nodes
ETAPSUM=zeros(Ny1,Nx1);
GGGPSUM=zeros(Ny1,Nx1);
SXXSUM=zeros(Ny1,Nx1);
RHOSUM=zeros(Ny1,Nx1);
WTPSUM=zeros(Ny1,Nx1);

for m=1:1:marknum
% Air/Lithosphere/channel
if(ym(m)<=ysize-xm(m)*tan((aplate+atopo)/180*pi))
% Air
if(tm(m)==2)
tm(m)=3;
end
elseif(ym(m)>yleft+(xm(m)-xleft)/(xright-xleft)*(yright-yleft) &&...
xm(m) % Channel
tm(m)=2;
else
if(ym(m)>ysize-xm(m)*tan(aplate/180*pi))
% Lithosphere
tm(m)=1;
end
end

% Interpolation to basic nodes
% Define i,j indexes for the upper left node
j=fix((xm(m)-x(1))/dx)+1;
i=fix((ym(m)-y(1))/dy)+1;
if(j<1)
j=1;
elseif(j>Nx-1)
j=Nx-1;
end
if(i<1)
i=1;
elseif(i>Ny-1)
i=Ny-1;
end
% Compute distances
dxmj=xm(m)-x(j);
dymi=ym(m)-y(i);
% Compute weights
wtmij=(1-dxmj/dx)*(1-dymi/dy);
wtmi1j=(1-dxmj/dx)*(dymi/dy);
wtmij1=(dxmj/dx)*(1-dymi/dy);
wtmi1j1=(dxmj/dx)*(dymi/dy);
% Update properties
% i,j Node
ETASUM(i,j)=ETASUM(i,j)+etam(tm(m))*wtmij;
GGGSUM(i,j)=GGGSUM(i,j)+1/gggm(tm(m))*wtmij;
SXYSUM(i,j)=SXYSUM(i,j)+sxym(m)*wtmij;
WTSUM(i,j)=WTSUM(i,j)+wtmij;
% i+1,j Node
ETASUM(i+1,j)=ETASUM(i+1,j)+etam(tm(m))*wtmi1j;
GGGSUM(i+1,j)=GGGSUM(i+1,j)+1/gggm(tm(m))*wtmi1j;
SXYSUM(i+1,j)=SXYSUM(i+1,j)+sxym(m)*wtmi1j;
WTSUM(i+1,j)=WTSUM(i+1,j)+wtmi1j;
% i,j+1 Node
ETASUM(i,j+1)=ETASUM(i,j+1)+etam(tm(m))*wtmij1;
GGGSUM(i,j+1)=GGGSUM(i,j+1)+1/gggm(tm(m))*wtmij1;
SXYSUM(i,j+1)=SXYSUM(i,j+1)+sxym(m)*wtmij1;
WTSUM(i,j+1)=WTSUM(i,j+1)+wtmij1;
% i+1,j+1 Node
ETASUM(i+1,j+1)=ETASUM(i+1,j+1)+etam(tm(m))*wtmi1j1;
GGGSUM(i+1,j+1)=GGGSUM(i+1,j+1)+1/gggm(tm(m))*wtmi1j1;
SXYSUM(i+1,j+1)=SXYSUM(i+1,j+1)+sxym(m)*wtmi1j1;
WTSUM(i+1,j+1)=WTSUM(i+1,j+1)+wtmi1j1;

% Interpolation to vx-nodes
% Define i,j indexes for the upper left node
j=fix((xm(m)-xvx(1))/dx)+1;
i=fix((ym(m)-yvx(1))/dy)+1;
if(j<1)
j=1;
elseif(j>Nx-1)
j=Nx-1;
end
if(i<1)
i=1;
elseif(i>Ny)
i=Ny;
end
% Compute distances
dxmj=xm(m)-xvx(j);
dymi=ym(m)-yvx(i);
% Compute weights
wtmij=(1-dxmj/dx)*(1-dymi/dy);
wtmi1j=(1-dxmj/dx)*(dymi/dy);
wtmij1=(dxmj/dx)*(1-dymi/dy);
wtmi1j1=(dxmj/dx)*(dymi/dy);
% Update properties
% i,j Node
RHOXSUM(i,j)=RHOXSUM(i,j)+rhom(tm(m))*wtmij;
WTXSUM(i,j)=WTXSUM(i,j)+wtmij;
% i+1,j Node
RHOXSUM(i+1,j)=RHOXSUM(i+1,j)+rhom(tm(m))*wtmi1j;
WTXSUM(i+1,j)=WTXSUM(i+1,j)+wtmi1j;
% i,j+1 Node
RHOXSUM(i,j+1)=RHOXSUM(i,j+1)+rhom(tm(m))*wtmij1;
WTXSUM(i,j+1)=WTXSUM(i,j+1)+wtmij1;
% i+1,j+1 Node
RHOXSUM(i+1,j+1)=RHOXSUM(i+1,j+1)+rhom(tm(m))*wtmi1j1;
WTXSUM(i+1,j+1)=WTXSUM(i+1,j+1)+wtmi1j1;

% Interpolation to vy-nodes
% Define i,j indexes for the upper left node
j=fix((xm(m)-xvy(1))/dx)+1;
i=fix((ym(m)-yvy(1))/dy)+1;
if(j<1)
j=1;
elseif(j>Nx)
j=Nx;
end
if(i<1)
i=1;
elseif(i>Ny-1)
i=Ny-1;
end
% Compute distances
dxmj=xm(m)-xvy(j);
dymi=ym(m)-yvy(i);
% Compute weights
wtmij=(1-dxmj/dx)*(1-dymi/dy);
wtmi1j=(1-dxmj/dx)*(dymi/dy);
wtmij1=(dxmj/dx)*(1-dymi/dy);
wtmi1j1=(dxmj/dx)*(dymi/dy);
% Update properties
% i,j Node
RHOYSUM(i,j)=RHOYSUM(i,j)+rhom(tm(m))*wtmij;
WTYSUM(i,j)=WTYSUM(i,j)+wtmij;
% i+1,j Node
RHOYSUM(i+1,j)=RHOYSUM(i+1,j)+rhom(tm(m))*wtmi1j;
WTYSUM(i+1,j)=WTYSUM(i+1,j)+wtmi1j;
% i,j+1 Node
RHOYSUM(i,j+1)=RHOYSUM(i,j+1)+rhom(tm(m))*wtmij1;
WTYSUM(i,j+1)=WTYSUM(i,j+1)+wtmij1;
% i+1,j+1 Node
RHOYSUM(i+1,j+1)=RHOYSUM(i+1,j+1)+rhom(tm(m))*wtmi1j1;
WTYSUM(i+1,j+1)=WTYSUM(i+1,j+1)+wtmi1j1;

% Interpolation to P-nodes
% Define i,j indexes for the upper left node
j=fix((xm(m)-xp(1))/dx)+1;
i=fix((ym(m)-yp(1))/dy)+1;
if(j<1)
j=1;
elseif(j>Nx)
j=Nx;
end
if(i<1)
i=1;
elseif(i>Ny)
i=Ny;
end
% Compute distances
dxmj=xm(m)-xp(j);
dymi=ym(m)-yp(i);
% Compute weights
wtmij=(1-dxmj/dx)*(1-dymi/dy);
wtmi1j=(1-dxmj/dx)*(dymi/dy);
wtmij1=(dxmj/dx)*(1-dymi/dy);
wtmi1j1=(dxmj/dx)*(dymi/dy);
% Update properties
% i,j Node
ETAPSUM(i,j)=ETAPSUM(i,j)+etam(tm(m))*wtmij;
GGGPSUM(i,j)=GGGPSUM(i,j)+1/gggm(tm(m))*wtmij;
SXXSUM(i,j)=SXXSUM(i,j)+sxxm(m)*wtmij;
RHOSUM(i,j)=RHOSUM(i,j)+rhom(tm(m))*wtmij;
WTPSUM(i,j)=WTPSUM(i,j)+wtmij;
% i+1,j Node
ETAPSUM(i+1,j)=ETAPSUM(i+1,j)+etam(tm(m))*wtmi1j;
GGGPSUM(i+1,j)=GGGPSUM(i+1,j)+1/gggm(tm(m))*wtmi1j;
SXXSUM(i+1,j)=SXXSUM(i+1,j)+sxxm(m)*wtmi1j;
RHOSUM(i+1,j)=RHOSUM(i+1,j)+rhom(tm(m))*wtmi1j;
WTPSUM(i+1,j)=WTPSUM(i+1,j)+wtmi1j;
% i,j+1 Node
ETAPSUM(i,j+1)=ETAPSUM(i,j+1)+etam(tm(m))*wtmij1;
GGGPSUM(i,j+1)=GGGPSUM(i,j+1)+1/gggm(tm(m))*wtmij1;
SXXSUM(i,j+1)=SXXSUM(i,j+1)+sxxm(m)*wtmij1;
RHOSUM(i,j+1)=RHOSUM(i,j+1)+rhom(tm(m))*wtmij1;
WTPSUM(i,j+1)=WTPSUM(i,j+1)+wtmij1;
% i+1,j+1 Node
ETAPSUM(i+1,j+1)=ETAPSUM(i+1,j+1)+etam(tm(m))*wtmi1j1;
GGGPSUM(i+1,j+1)=GGGPSUM(i+1,j+1)+1/gggm(tm(m))*wtmi1j1;
SXXSUM(i+1,j+1)=SXXSUM(i+1,j+1)+sxxm(m)*wtmi1j1;
RHOSUM(i+1,j+1)=RHOSUM(i+1,j+1)+rhom(tm(m))*wtmi1j1;
WTPSUM(i+1,j+1)=WTPSUM(i+1,j+1)+wtmi1j1;
end
% Compute physical properties
% Basic nodes
for j=1:1:Nx
for i=1:1:Ny
if(WTSUM(i,j)>0)
ETA(i,j)=ETASUM(i,j)/WTSUM(i,j);
GGG(i,j)=1/(GGGSUM(i,j)/WTSUM(i,j));
SXY0(i,j)=SXYSUM(i,j)/WTSUM(i,j);
end
end
end
% Vx-nodes
for j=1:1:Nx1
for i=1:1:Ny1
if(WTXSUM(i,j)>0)
RHOX(i,j)=RHOXSUM(i,j)/WTXSUM(i,j);
end
end
end
% Vy-nodes
for j=1:1:Nx1
for i=1:1:Ny1
if(WTYSUM(i,j)>0)
RHOY(i,j)=RHOYSUM(i,j)/WTYSUM(i,j);
end
end
end
% P-nodes
for j=1:1:Nx1
for i=1:1:Ny1
if(WTPSUM(i,j)>0)
ETAP(i,j)=ETAPSUM(i,j)/WTPSUM(i,j);
GGGP(i,j)=1/(GGGPSUM(i,j)/WTPSUM(i,j));
SXX0(i,j)=SXXSUM(i,j)/WTPSUM(i,j);
RHO(i,j)=RHOSUM(i,j)/WTPSUM(i,j);
end
end
end

% Mechanical Solution
% Composing global matrixes L(), R() for Stokes and continuity equations
for j=1:1:Nx1
for i=1:1:Ny1
% Define global indexes in algebraic space
kvx=((j-1)*Ny1+i-1)*3+1; % Vx
kvy=kvx+1; % Vy
kpm=kvx+2; % P

% Vx equation External points
if(i==1 || i==Ny1 || j==1 || j==Nx || j==Nx1)
% Boundary Condition
% 1*Vx=0
L(kvx,kvx)=1; % Left part
R(kvx)=0; % Right part
% Top boundary
if(i==1 && j>1 && j L(kvx,kvx+3)=bctop; % Left part
end
% Bottom boundary
if(i==Ny1 && j>1 && j L(kvx,kvx-3)=bcbottom; % Left part
% Plate velocity
if(j R(kvx)=vplate*2; % Right part
end
end
else
% Internal points: x-Stokes eq.
% Vx2
% |
% Vy1 | Vy3
% |
% Vx1-P1-Vx3-P2-Vx5
% |
% Vy2 | Vy4
% |
% Vx4
%
% Computational viscosity
ETA1=ETA(i-1,j)*GGG(i-1,j)*dt/(GGG(i-1,j)*dt+ETA(i-1,j));
ETA2=ETA(i,j)*GGG(i,j)*dt/(GGG(i,j)*dt+ETA(i,j));
ETAP1=ETAP(i,j)*GGGP(i,j)*dt/(GGGP(i,j)*dt+ETAP(i,j));
ETAP2=ETAP(i,j+1)*GGGP(i,j+1)*dt/(GGGP(i,j+1)*dt+ETAP(i,j+1));
% Old stresses
SXY1=SXY0(i-1,j)*ETA(i-1,j)/(GGG(i-1,j)*dt+ETA(i-1,j));
SXY2=SXY0(i,j)*ETA(i,j)/(GGG(i,j)*dt+ETA(i,j));
SXX1=SXX0(i,j)*ETAP(i,j)/(GGGP(i,j)*dt+ETAP(i,j));
SXX2=SXX0(i,j+1)*ETAP(i,j+1)/(GGGP(i,j+1)*dt+ETAP(i,j+1));
% Density gradients
dRHOdx=(RHOX(i,j+1)-RHOX(i,j-1))/2/dx;
dRHOdy=(RHOX(i+1,j)-RHOX(i-1,j))/2/dy;
% Left part
L(kvx,kvx-Ny1*3)=2*ETAP1/dx^2; % Vx1
L(kvx,kvx-3)=ETA1/dy^2; % Vx2
L(kvx,kvx)=-2*(ETAP1+ETAP2)/dx^2-...
(ETA1+ETA2)/dy^2-...
dRHOdx*gx*dt; % Vx3
L(kvx,kvx+3)=ETA2/dy^2; % Vx4
L(kvx,kvx+Ny1*3)=2*ETAP2/dx^2; % Vx5
L(kvx,kvy)=-ETA2/dx/dy-dRHOdy*gx*dt/4; % Vy2
L(kvx,kvy+Ny1*3)=ETA2/dx/dy-dRHOdy*gx*dt/4; % Vy4
L(kvx,kvy-3)=ETA1/dx/dy-dRHOdy*gx*dt/4; % Vy1
L(kvx,kvy+Ny1*3-3)=-ETA1/dx/dy-dRHOdy*gx*dt/4; % Vy3
L(kvx,kpm)=pscale/dx; % P1
L(kvx,kpm+Ny1*3)=-pscale/dx; % P2
% Right part
R(kvx)=-RHOX(i,j)*gx-(SXY2-SXY1)/dy-(SXX2-SXX1)/dx;
end

% Vy equation External points
if(j==1 || j==Nx1 || i==1 || i==Ny || i==Ny1)
% Boundary Condition
% 1*Vy=0
L(kvy,kvy)=1; % Left part
R(kvy)=0; % Right part
% Left boundary
if(j==1 && i>1 && i L(kvy,kvy+3*Ny1)=bcleft; % Left part
end
% Right boundary
if(j==Nx1 && i>1 && i L(kvy,kvy-3*Ny1)=bcright; % Left part
end
else
% Internal points: y-Stokes eq.
% Vy2
% |
% Vx1 P1 Vx3
% |
% Vy1----Vy3----Vy5
% |
% Vx2 P2 Vx4
% |
% Vy4
%
% Computational viscosity
ETA1=ETA(i,j-1)*GGG(i,j-1)*dt/(GGG(i,j-1)*dt+ETA(i,j-1));
ETA2=ETA(i,j)*GGG(i,j)*dt/(GGG(i,j)*dt+ETA(i,j));
ETAP1=ETAP(i,j)*GGGP(i,j)*dt/(GGGP(i,j)*dt+ETAP(i,j));
ETAP2=ETAP(i+1,j)*GGGP(i+1,j)*dt/(GGGP(i+1,j)*dt+ETAP(i+1,j));
% Old stresses
SXY1=SXY0(i,j-1)*ETA(i,j-1)/(GGG(i,j-1)*dt+ETA(i,j-1));
SXY2=SXY0(i,j)*ETA(i,j)/(GGG(i,j)*dt+ETA(i,j));
SYY1=-SXX0(i,j)*ETAP(i,j)/(GGGP(i,j)*dt+ETAP(i,j));
SYY2=-SXX0(i+1,j)*ETAP(i+1,j)/(GGGP(i+1,j)*dt+ETAP(i+1,j));
% Density gradients
dRHOdx=(RHOY(i,j+1)-RHOY(i,j-1))/2/dx;
dRHOdy=(RHOY(i+1,j)-RHOY(i-1,j))/2/dy;
% Left part
L(kvy,kvy-Ny1*3)=ETA1/dx^2; % Vy1
L(kvy,kvy-3)=2*ETAP1/dy^2; % Vy2
L(kvy,kvy)=-2*(ETAP1+ETAP2)/dy^2-...
(ETA1+ETA2)/dx^2-...
dRHOdy*gy*dt; % Vy3
L(kvy,kvy+3)=2*ETAP2/dy^2; % Vy4
L(kvy,kvy+Ny1*3)=ETA2/dx^2; % Vy5
L(kvy,kvx)=-ETA2/dx/dy-dRHOdx*gy*dt/4; %Vx3
L(kvy,kvx+3)=ETA2/dx/dy-dRHOdx*gy*dt/4; %Vx4
L(kvy,kvx-Ny1*3)=ETA1/dx/dy-dRHOdx*gy*dt/4; %Vx1
L(kvy,kvx+3-Ny1*3)=-ETA1/dx/dy-dRHOdx*gy*dt/4; %Vx2
L(kvy,kpm)=pscale/dy; % P1
L(kvy,kpm+3)=-pscale/dy; % P2

% Right part
R(kvy)=-RHOY(i,j)*gy-(SXY2-SXY1)/dx-(SYY2-SYY1)/dy;
end

% P equation External points
if(i==1 || j==1 || i==Ny1 || j==Nx1 ||...
(i==2 && j==2))
% Boundary Condition
% 1*P=0
L(kpm,kpm)=1; % Left part
R(kpm)=0; % Right part
% Real BC
if(i==2 && j==2)
L(kpm,kpm)=1*pscale; %Left part
R(kpm)=1e+5; % Right part
end
else
% Internal points: continuity eq.
% dVx/dx+dVy/dy=0
% Vy1
% |
% Vx1--P--Vx2
% |
% Vy2
%
% Left part
L(kpm,kvx-Ny1*3)=-1/dx; % Vx1
L(kpm,kvx)=1/dx; % Vx2
L(kpm,kvy-3)=-1/dy; % Vy1
L(kpm,kvy)=1/dy; % Vy2
% Right part
R(kpm)=0;
end

end
end

% 4) Solving matrixes, reloading solution
S=L\R; % Obtaining algebraic vector of solutions S()
% Reload solutions S() to vx(), vy(), p()
% Going through all grid points
for j=1:1:Nx1
for i=1:1:Ny1
% Define global indexes in algebraic space
kvx=((j-1)*Ny1+i-1)*3+1; % Vx
kvy=kvx+1; % Vy
kpm=kvx+2; % P
% Reload solution
vx(i,j)=S(kvx);
vy(i,j)=S(kvy);
pr(i,j)=S(kpm)*pscale;
end
end

% Define displacement timestep dtm
dtm=dt;
maxvx=max(max(abs(vx)));
maxvy=max(max(abs(vy)));
if(dtm*maxvx>dxymax*dx)
dtm=dxymax*dx/maxvx;
end
if(dtm*maxvy>dxymax*dy)
dtm=dxymax*dy/maxvy;
end

% Compute Stress, stress change and strain rate components
% Compute EPSILONxy, SIGMAxy in basic nodes
EXY=zeros(Ny,Nx); % Strain rate EPSILONxy, 1/s
SXY=zeros(Ny,Nx); % Stress SIGMAxy, Pa
DSXY=zeros(Ny,Nx); % Stress change SIGMAxy, Pa
for j=1:1:Nx
for i=1:1:Ny
% EXY,SXY, DSXY
EXY(i,j)=0.5*((vx(i+1,j)-vx(i,j))/dy+...
(vy(i,j+1)-vy(i,j))/dx);
SXY(i,j)=2*ETA(i,j)*EXY(i,j)*GGG(i,j)*dtm/(GGG(i,j)*dtm+ETA(i,j))+...
SXY0(i,j)*ETA(i,j)/(GGG(i,j)*dtm+ETA(i,j));
DSXY(i,j)=SXY(i,j)-SXY0(i,j);
end
end
% Compute EPSILONxx, SIGMA'xx in pressure nodes
EXX=zeros(Ny1,Nx1); % Strain rate EPSILONxx, 1/s
SXX=zeros(Ny1,Nx1); % Stress SIGMA'xx, Pa
DSXX=zeros(Ny1,Nx1); % Stress change SIGMA'xx, Pa
for j=2:1:Nx
for i=2:1:Ny
% EXX
EXX(i,j)=(vx(i,j)-vx(i,j-1))/dx;
% SXX
SXX(i,j)=2*ETAP(i,j)*EXX(i,j)*GGGP(i,j)*dtm/(GGGP(i,j)*dtm+ETAP(i,j))+...
SXX0(i,j)*ETAP(i,j)/(GGGP(i,j)*dtm+ETAP(i,j));
DSXX(i,j)=SXX(i,j)-SXX0(i,j);
end
end

% Save nodal stress changes
DSXX0=DSXX;
DSXY0=DSXY;
% Apply subgrid stress diffusion on markers
if(dsubgrids>0)
SXYSUM=zeros(Ny,Nx);
WTSUM=zeros(Ny,Nx);
SXXSUM=zeros(Ny1,Nx1);
WTPSUM=zeros(Ny1,Nx1);
for m=1:1:marknum
% SIGMA'xx
% Define i,j indexes for the upper left node
j=fix((xm(m)-xp(1))/dx)+1;
i=fix((ym(m)-yp(1))/dy)+1;
if(j<1)
j=1;
elseif(j>Nx)
j=Nx;
end
if(i<1)
i=1;
elseif(i>Ny)
i=Ny;
end
% Compute distances
dxmj=xm(m)-xp(j);
dymi=ym(m)-yp(i);
% Compute weights
wtmij=(1-dxmj/dx)*(1-dymi/dy);
wtmi1j=(1-dxmj/dx)*(dymi/dy);
wtmij1=(dxmj/dx)*(1-dymi/dy);
wtmi1j1=(dxmj/dx)*(dymi/dy);
% Compute marker-node SIGMA'xx difference
dsxxm0=sxxm(m)-(SXX0(i,j)*wtmij+SXX0(i+1,j)*wtmi1j+...
SXX0(i,j+1)*wtmij1+SXX0(i+1,j+1)*wtmi1j1);
% Relax stress difference
dsxxm1=dsxxm0*exp(-dsubgrids*dtm/(etam(tm(m))/gggm(tm(m))));
% Correct marker stress
ddsxxm=dsxxm1-dsxxm0;
sxxm(m)=sxxm(m)+ddsxxm;
% Update subgrid diffusion on nodes
% i,j Node
SXXSUM(i,j)=SXXSUM(i,j)+ddsxxm*wtmij;
WTPSUM(i,j)=WTPSUM(i,j)+wtmij;
% i+1,j Node
SXXSUM(i+1,j)=SXXSUM(i+1,j)+ddsxxm*wtmi1j;
WTPSUM(i+1,j)=WTPSUM(i+1,j)+wtmi1j;
% i,j+1 Node
SXXSUM(i,j+1)=SXXSUM(i,j+1)+ddsxxm*wtmij1;
WTPSUM(i,j+1)=WTPSUM(i,j+1)+wtmij1;
% i+1,j+1 Node
SXXSUM(i+1,j+1)=SXXSUM(i+1,j+1)+ddsxxm*wtmi1j1;
WTPSUM(i+1,j+1)=WTPSUM(i+1,j+1)+wtmi1j1;

% SIGMAxy
% Define i,j indexes for the upper left node
j=fix((xm(m)-x(1))/dx)+1;
i=fix((ym(m)-y(1))/dy)+1;
if(j<1)
j=1;
elseif(j>Nx-1)
j=Nx-1;
end
if(i<1)
i=1;
elseif(i>Ny-1)
i=Ny-1;
end
% Compute distances
dxmj=xm(m)-x(j);
dymi=ym(m)-y(i);
% Compute weights
wtmij=(1-dxmj/dx)*(1-dymi/dy);
wtmi1j=(1-dxmj/dx)*(dymi/dy);
wtmij1=(dxmj/dx)*(1-dymi/dy);
wtmi1j1=(dxmj/dx)*(dymi/dy);
% Compute marker-node SIGMAxy difference
dsxym0=sxym(m)-(SXY0(i,j)*wtmij+SXY0(i+1,j)*wtmi1j+...
SXY0(i,j+1)*wtmij1+SXY0(i+1,j+1)*wtmi1j1);
% Relax stress difference
dsxym1=dsxym0*exp(-dsubgrids*dtm/(etam(tm(m))/gggm(tm(m))));
% Correct marker stress
ddsxym=dsxym1-dsxym0;
sxym(m)=sxym(m)+ddsxym;
% Update subgrid diffusion on nodes
% i,j Node
SXYSUM(i,j)=SXYSUM(i,j)+ddsxym*wtmij;
WTSUM(i,j)=WTSUM(i,j)+wtmij;
% i+1,j Node
SXYSUM(i+1,j)=SXYSUM(i+1,j)+ddsxym*wtmi1j;
WTSUM(i+1,j)=WTSUM(i+1,j)+wtmi1j;
% i,j+1 Node
SXYSUM(i,j+1)=SXYSUM(i,j+1)+ddsxym*wtmij1;
WTSUM(i,j+1)=WTSUM(i,j+1)+wtmij1;
% i+1,j+1 Node
SXYSUM(i+1,j+1)=SXYSUM(i+1,j+1)+ddsxym*wtmi1j1;
WTSUM(i+1,j+1)=WTSUM(i+1,j+1)+wtmi1j1;
end
% Compute DSXXsubgrid
DSXXsubgrid=zeros(Ny1,Nx1);
% P-nodes
for j=2:1:Nx
for i=2:1:Ny
if(WTPSUM(i,j)>0)
DSXXsubgrid(i,j)=SXXSUM(i,j)/WTPSUM(i,j);
end
end
end
% Correct DSXX
DSXX=DSXX-DSXXsubgrid;
% Compute DSXYsubgrid
DSXYsubgrid=zeros(Ny,Nx);
% Basic nodes
for j=1:1:Nx
for i=1:1:Ny
if(WTSUM(i,j)>0)
DSXYsubgrid(i,j)=SXYSUM(i,j)/WTSUM(i,j);
end
end
end
% Correct DSXY
DSXY=DSXY-DSXYsubgrid;
end

% Interpolate DSXX, DSXY to markers
for m=1:1:marknum
% SIGMA'xx
% Define i,j indexes for the upper left node
j=fix((xm(m)-xp(1))/dx)+1;
i=fix((ym(m)-yp(1))/dy)+1;
if(j<2)
j=2;
elseif(j>Nx-1)
j=Nx-1;
end
if(i<2)
i=2;
elseif(i>Ny-1)
i=Ny-1;
end
% Compute distances
dxmj=xm(m)-xp(j);
dymi=ym(m)-yp(i);
% Compute weights
wtmij=(1-dxmj/dx)*(1-dymi/dy);
wtmi1j=(1-dxmj/dx)*(dymi/dy);
wtmij1=(dxmj/dx)*(1-dymi/dy);
wtmi1j1=(dxmj/dx)*(dymi/dy);
% Update marker by SIGMA'xx change
sxxm(m)=sxxm(m)+(DSXX(i,j)*wtmij+DSXX(i+1,j)*wtmi1j+...
DSXX(i,j+1)*wtmij1+DSXX(i+1,j+1)*wtmi1j1);

% SIGMAxy
% Define i,j indexes for the upper left node
j=fix((xm(m)-x(1))/dx)+1;
i=fix((ym(m)-y(1))/dy)+1;
if(j<1)
j=1;
elseif(j>Nx-1)
j=Nx-1;
end
if(i<1)
i=1;
elseif(i>Ny-1)
i=Ny-1;
end
% Compute distances
dxmj=xm(m)-x(j);
dymi=ym(m)-y(i);
% Compute weights
wtmij=(1-dxmj/dx)*(1-dymi/dy);
wtmi1j=(1-dxmj/dx)*(dymi/dy);
wtmij1=(dxmj/dx)*(1-dymi/dy);
wtmi1j1=(dxmj/dx)*(dymi/dy);
% Update marker by SIGMA'xx change
sxym(m)=sxym(m)+(DSXY(i,j)*wtmij+DSXY(i+1,j)*wtmi1j+...
DSXY(i,j+1)*wtmij1+DSXY(i+1,j+1)*wtmi1j1);
end

% Compute velocity in pressure nodes
% vx
for j=2:1:Nx
for i=2:1:Ny
vxp(i,j)=(vx(i,j)+vx(i,j-1))/2;
end
end
% Apply BC
% Top
vxp(1,2:Nx-1)=-bctop*vxp(2,2:Nx-1);
% Bottom
vxp(Ny1,2:Nx-1)=-bcbottom*vxp(Ny,2:Nx-1);
% Left
vxp(:,1)=-vxp(:,2);
% Right
vxp(:,Nx1)=-vxp(:,Nx);
% vy
for j=2:1:Nx
for i=2:1:Ny
vyp(i,j)=(vy(i,j)+vy(i-1,j))/2;
end
end
% Apply BC
% Left
vyp(2:Ny-1,1)=-bcleft*vyp(2:Ny-1,2);
% Right
vyp(2:Ny-1,Nx1)=-bcright*vyp(2:Ny-1,Nx); % Free slip
% Top
vyp(1,:)=-vyp(2,:);
% Bottom
vyp(Ny1,:)=-vyp(Ny,:);

% Compute rotation rate wyx=1/2(dVy/dx-dVx/dy) for basic nodes
for i=1:1:Ny
for j=1:1:Nx
wyx(i,j)=0.5*((vy(i,j+1)-vy(i,j))/dx-(vx(i+1,j)-vx(i,j))/dy);
end
end

% Move markers with 4th order Runge-Kutta
vxm=zeros(4,1);
vym=zeros(4,1);
for m=1:1:marknum
% Interpolate local rotation rate
% Define i,j indexes for the upper left node
j=fix((xm(m)-x(1))/dx)+1;
i=fix((ym(m)-y(1))/dy)+1;
if(j<1)
j=1;
elseif(j>Nx-1)
j=Nx-1;
end
if(i<1)
i=1;
elseif(i>Ny-1)
i=Ny-1;
end
% Compute distances
dxmj=xm(m)-x(j);
dymi=ym(m)-y(i);
% Compute weights
wtmij=(1-dxmj/dx)*(1-dymi/dy);
wtmi1j=(1-dxmj/dx)*(dymi/dy);
wtmij1=(dxmj/dx)*(1-dymi/dy);
wtmi1j1=(dxmj/dx)*(dymi/dy);
% Compute vx velocity
omegam=wyx(i,j)*wtmij+wyx(i+1,j)*wtmi1j+...
wyx(i,j+1)*wtmij1+wyx(i+1,j+1)*wtmi1j1;
% Analytical stress rotation using SIGMA'xx=-SIGMA'yy
THETA=dtm*omegam; % Incremental rotation angle
sxxmnew=sxxm(m)*cos(THETA)^2-sxxm(m)*sin(THETA)^2-sxym(m)*sin(2*THETA);
sxymnew=sxxm(m)*sin(2*THETA)+sxym(m)*cos(2*THETA);
sxxm(m)=sxxmnew; sxym(m)=sxymnew;

% Save initial marker's coordinates
xA=xm(m);
yA=ym(m);
for rk=1:1:4
% Interpolate vxp,vyp
% Define i,j indexes for the upper left node
j=fix((xm(m)-xp(1))/dx)+1;
i=fix((ym(m)-yp(1))/dy)+1;
if(j<1)
j=1;
elseif(j>Nx)
j=Nx;
end
if(i<1)
i=1;
elseif(i>Ny)
i=Ny;
end
% Compute distances
dxmj=xm(m)-xp(j);
dymi=ym(m)-yp(i);
% Compute weights
wtmij=(1-dxmj/dx)*(1-dymi/dy);
wtmi1j=(1-dxmj/dx)*(dymi/dy);
wtmij1=(dxmj/dx)*(1-dymi/dy);
wtmi1j1=(dxmj/dx)*(dymi/dy);
% Compute vx, vy velocity
vxm(rk)=vxp(i,j)*wtmij+vxp(i+1,j)*wtmi1j+...
vxp(i,j+1)*wtmij1+vxp(i+1,j+1)*wtmi1j1;
vym(rk)=vyp(i,j)*wtmij+vyp(i+1,j)*wtmi1j+...
vyp(i,j+1)*wtmij1+vyp(i+1,j+1)*wtmi1j1;

% Interpolate vx
% Define i,j indexes for the upper left node
j=fix((xm(m)-xvx(1))/dx)+1;
i=fix((ym(m)-yvx(1))/dy)+1;
if(j<1)
j=1;
elseif(j>Nx-1)
j=Nx-1;
end
if(i<1)
i=1;
elseif(i>Ny)
i=Ny;
end
% Compute distances
dxmj=xm(m)-xvx(j);
dymi=ym(m)-yvx(i);
% Compute weights
wtmij=(1-dxmj/dx)*(1-dymi/dy);
wtmi1j=(1-dxmj/dx)*(dymi/dy);
wtmij1=(dxmj/dx)*(1-dymi/dy);
wtmi1j1=(dxmj/dx)*(dymi/dy);
% Compute vx velocity
vxm(rk)=vpratio*vxm(rk)+(1-vpratio)*(vx(i,j)*wtmij+vx(i+1,j)*wtmi1j+...
vx(i,j+1)*wtmij1+vx(i+1,j+1)*wtmi1j1);

% Interpolate vy
% Define i,j indexes for the upper left node
j=fix((xm(m)-xvy(1))/dx)+1;
i=fix((ym(m)-yvy(1))/dy)+1;
if(j<1)
j=1;
elseif(j>Nx)
j=Nx;
end
if(i<1)
i=1;
elseif(i>Ny-1)
i=Ny-1;
end
% Compute distances
dxmj=xm(m)-xvy(j);
dymi=ym(m)-yvy(i);
% Compute weights
wtmij=(1-dxmj/dx)*(1-dymi/dy);
wtmi1j=(1-dxmj/dx)*(dymi/dy);
wtmij1=(dxmj/dx)*(1-dymi/dy);
wtmi1j1=(dxmj/dx)*(dymi/dy);
% Compute vx velocity
vym(rk)=vpratio*vym(rk)+(1-vpratio)*(vy(i,j)*wtmij+vy(i+1,j)*wtmi1j+...
vy(i,j+1)*wtmij1+vy(i+1,j+1)*wtmi1j1);

% Change coordinates to obtain B,C,D points
if(rk==1 || rk==2)
xm(m)=xA+dtm/2*vxm(rk);
ym(m)=yA+dtm/2*vym(rk);
elseif(rk==3)
xm(m)=xA+dtm*vxm(rk);
ym(m)=yA+dtm*vym(rk);
end
end
% Restore initial coordinates
xm(m)=xA;
ym(m)=yA;
% Compute effective velocity
vxmeff=1/6*(vxm(1)+2*vxm(2)+2*vxm(3)+vxm(4));
vymeff=1/6*(vym(1)+2*vym(2)+2*vym(3)+vym(4));
% Move markers
xm(m)=xm(m)+dtm*vxmeff;
ym(m)=ym(m)+dtm*vymeff;
end

% Update timesum
timesum=timesum+dtm;

figure(1);colormap('Jet');clf
subplot(4,2,1)
pcolor(x/1000,y/1000,log10(ETA));
shading flat;
axis ij image;
colorbar
title('log10ETA, Pa*s')
hold on
quiver(xp(3:5:Nx1)/1000,yp(3:5:Ny1)/1000,vxp(3:5:Ny,3:5:Nx1),vyp(3:5:Ny1,3:5:Nx1),'w')

subplot(4,2,2)
pcolor(xp/1000,yp/1000,RHO)
shading interp;
axis ij image;
colorbar
title('Density, kg/m^3')

subplot(4,2,3)
pcolor(xp/1000,yp/1000,vxp)
shading interp;
axis ij image;
colorbar
title(['Vx, m/s time=',num2str(timesum/(365.25*24*3600)),' Yr'])
hold on
quiver(xp(3:5:Nx1)/1000,yp(3:5:Ny1)/1000,vxp(3:5:Ny,3:5:Nx1),vyp(3:5:Ny1,3:5:Nx1),'k')

subplot(4,2,4)
pcolor(xp/1000,yp/1000,vyp)
shading interp;
axis ij image;
colorbar
title('Vy, m/s')
hold on
quiver(xp(3:5:Nx1)/1000,yp(3:5:Ny1)/1000,vxp(3:5:Ny,3:5:Nx1),vyp(3:5:Ny1,3:5:Nx1),'k')

subplot(4,2,5)
pcolor(xp(2:Nx)/1000,yp(2:Ny)/1000,pr(2:Ny,2:Nx))
shading interp;
axis ij image;
colorbar
title('Pressure, Pa')
hold on
quiver(xp(3:5:Nx1)/1000,yp(3:5:Ny1)/1000,vxp(3:5:Ny,3:5:Nx1),vyp(3:5:Ny1,3:5:Nx1),'k')

subplot(4,2,6)
pcolor(xp(2:Nx)/1000,yp(2:Ny)/1000,pr(2:Ny,2:Nx)-prl(2:Ny,2:Nx))
% pcolor(xp(2:Nx)/1000,yp(2:Ny)/1000,prl(2:Ny,2:Nx))
shading interp;
axis ij image;
colorbar
title('Overpressure, Pa')

subplot(4,2,7)
pcolor(xp/1000,yp/1000,SXX)
shading interp;
axis ij image;
colorbar
title('TAUxx, Pa')

subplot(4,2,8)
pcolor(x/1000,y/1000,SXY)
shading interp;
axis ij image;
colorbar
title('TAUxy, Pa')

% figure(2);colormap('Jet');clf
% plot(xm,ym,'. k');
% axis ij image

pause(0.1)

end

Hide

AR by John Platt on behalf of the Authors (24 Dec 2018) Author's response Manuscript

ED: Referee Nomination & Report Request started (03 Jan 2019) by Taras Gerya

RR by Stefan Markus Schmalholz (28 Jan 2019)

Suggestions for revision or reasons for rejection

I find the revised version of the comment by Platt now acceptable for publication. After the online discussions (which are in my opinion useful and interesting) related to the first version of the comment of Platt I am looking forward to the reply of Marques and co-authors. I have only a few comments to the revised comment of Platt.

Line 22-23, page 3: I find this sentence misleading. If the resistance of the plate becomes larger, then the corresponding deflection, due to a given normal load, becomes smaller. Therefore, an extremely small deflection can correspond to an extremely large resistance. The boundary condition of a straight channel wall with zero deflection applied by Marques et al. (2018) can be justified by two assumptions: (i) the deflection is small with respect to the length of the channel wall so that it is negligible. (ii) the straight channel wall corresponds to the elastically-deflected state implying that without dynamic pressure the upper channel wall would be curved.

Line 2-11, page 5: This is a selective view on the topic of the strength of the lithosphere. There are many studies suggesting much higher shear stresses in a deforming continental lithosphere, especially when flexure is considered. Many studies on the strength of the lithosphere considering so-called yield-strength envelopes, based on flow laws from rock deformation experiments, exhibit one or more levels in the lithosphere with stresses corresponding to shear stresses significantly larger than 120 MPa (e.g. Burov, 2011, and references therein). For example, if strength in the hanging wall in depths larger than 40 km would be dominated by diabase then the stress could be significantly larger than 120 MPa (e.g. Burov, 2011, his figure 7). Furthermore, flow laws for dry anorthite indicate that stresses in the lower crust could be significantly larger than 120 MPa (e.g. Rybacki et al., 2006).

Line 14-16, page 5: This suggestion that “the values of dynamic overpressure calculated by M2018 are at least an order of magnitude too high” is based on the assumption that the rocks forming the hanging wall exhibit shear stresses smaller than ca. 120 MPa. Marques et al. (2018) are aware of the fact that the magnitude of the dynamic pressure depends on the effective strength of the hanging wall, and they write in their manuscript: “Ultimately, tectonic overpressure (TOP) can only exist if the channel walls are strong enough.”. The effective strength of the hanging wall is a key parameter controlling the magnitude of dynamic pressure in the model of Marques et al. (2018). I agree that a dynamic pressure of 1.5 GPa within large channel regions, as suggested by Marques et al. (2018), is very difficult to justify. However, a magnitude of the dynamic pressure of 0.3 to 0.5 GPa is a factor of 5 to 3 times smaller than the value of 1.5 GPa and I think that values between 0.3 and 0.5 GPa for dynamic pressure cannot be excluded based on currently available flow laws for the lower crust. Hence, I still think that the above “suggestion” that values of dynamic overpressure are “at least an order of magnitude too high” is also very difficult to justify.

Best regards,
Stefan Schmalholz

References
Burov, E.B., 2011. Rheology and strength of the lithosphere. Marine and Petroleum Geology 28, 1402-1443.
Rybacki, E., Gottschalk, M., Wirth, R., Dresen, G., 2006. Influence of water fugacity and activation volume on the flow properties of fine-grained anorthite aggregates. Journal of Geophysical Research-Solid Earth 111.

Hide

RR by Evangelos Moulas (04 Feb 2019)

ED: Publish subject to minor revisions (review by editor) (05 Feb 2019) by Taras Gerya

AR by John Platt on behalf of the Authors (10 Feb 2019) Author's response Manuscript

ED: Publish as is (18 Feb 2019) by Taras Gerya

ED: Publish as is (18 Feb 2019) by Federico Rossetti (Executive editor)

AR by John Platt on behalf of the Authors (19 Feb 2019) Manuscript

Short summary

The channel flow model proposed by Marques et al (2018) for the Himalayas has a geometry that would not generate any excess pressure. The excess pressure calculated by the authors, based on a different and highly improbable geometry, is so high that the overlying rocks would not be able to contain it: they would bend or break in such a way as to relieve the pressure. The conclusions drawn by the authors are therefore unwarranted and misleading.