Dear glaciologists,

I have been working on estimating glacier thickness along flowline with only input information being mass balance, width, and surface topography. I assume steady state. The idea is similar to Farinotti et al.(2009), J. of Glac., but I don't want 2-D thickness, I just want it along flowline. Also, I assume I know mass balance along each point on the flowline and I am not calibrating to a known glacier volume.

I have found a paper written by Russian mathematicians Salamatin and Mazo in 1985 that has done what I want to do, theoretically. However, I don't speak Russian so that is a problem. I have attached the paper, but I don't think I can attach it to the list so if you'd like to see it, please email me at

[hidden email].

The problem I am having is this: I would like to solve the following fully implicit ODE for u(x)

0=k*(u(x)-z(x))* F1(s)+((u(x)-z(x))^2)* F2(s) -A(x) from 0<x<1

s=(u(x)-z(x))*abs(du/dx) , and F1(s)=s^((n+1)/2), F2(s)=(s^n)/(n+2)

with boundary condition u(1)=z(1)

This describes 0=sliding+ creep - all mass balance flux above.

We say we know z(x), A(x), and n, and we would like to solve for u(x). I cannot get this equation to not have singular iteration matrices when I use DAE solvers in fortran. I have tried using dassl.f (the latest one) using BDF methods with direct and preconditioned Krylov linear solvers; and coldae.f (the newest of colsys and colnew) using collocation, projected on constraint manifold. I got them off

www.netlib.org/ode/
The code will work in matlab with ode15i, but I'm thinking now that is because that code made more approximations than the fortran codes, and so it wasn't going singular at the boundary point x=1. Even so, the ode15i.m code crashed quite a bit when trying to get an initial point on some glaciers.

I'm wondering if anyone has any advice on a different code to solve this ODE, or, if you can read Russian, can you tell me if the section from page 100 to 101 tells me how the original authors solved the problem? In another paper translated into English, they say "the boundary-value problem . . . can be solved numerically without difficulty (reference to the paper only in Russian."

Thank you for your help,

Ashley Van Beuskom

_______________________________________________

You're subscribed to the CRYOLIST mailing list

To change your subscription options, visit

http://cryolist.org/To send a message to the list, email

[hidden email]