function [flux] = findMeltwaterFlux(porosity,alpha,kappa,effpressure,z,delta)

dNdz = gradient(effpressure)./gradient(z);

flux = (kappa.*porosity.^alpha).*(-1+delta.*dNdz);

end