function [ U ] = isostaticmassflux1D(x,t,varargin) % [ U ] = isostaticmassflux1D(x,t) % [ U ] = isostaticmassflux1D(x,t,water) % % Numerically determine the mass flux in the asthenosphere for an across %axis profile deriving from isostasy driven flow with differential thermal %subsidence % % x = location vector (assume km, but units are actually arbirtrary) % t = thermal age of lithosphere (governs subsidence rate) % water = optional binary vector defining which of locations are % submarine (1) and which are subaerial (0). Default is all submarine. % % U = mass flux (km^3/Myr/km of ridge length) % % Generalization of the mathematics in Conder (2012), Lithosphere % Developed by James Conder, March 2021 % % Included as supplemental materials to % Conder, J.A., Oceanic isostasy as a trigger to the rift to drift % transition, submitted to Geology, 2021 nx = length(x); cstar = 0.11; % km/sqrt(Myr) t = max(0.02,t); massgain = cstar./(2*sqrt(t)); %mass gain per area (km/Myr) if nargin > 2 submarine = varargin{1}; if length(submarine) == nx massgain = massgain.*submarine; %no ocean = no mass gain end end dx = (x(3:end) - x(1:end-2))/2; dx = [(x(2)-x(1)) dx (x(end)-x(end-1))]; meanmassgain = sum(dx.*massgain)/(x(end)-x(1)); dUdx = massgain - meanmassgain; U = 0*x; U(1) = dUdx(1)*dx(1); for i = 2:length(x) U(i) = U(i-1) + dUdx(i)*dx(i); end end