I'm trying to recreate the graphs in Figure 6.2 from this source (page 88). They get these plots by using the following solutions for velocity and temperature (both obtained from solving the basic parallel sided slab on slope problem).

$v_x = \frac{2 A (\rho g sin(\alpha))^{n}}{n+1}(H^{n+1} -(H-z)^{n+1})$
$T = T_s + \frac{q^{\perp}_{geo}}{\kappa}(H-z)+\frac{2AH^{n+3} (\rho g sin(\alpha))^{n+1}}{\kappa (n+2)} [1-\frac{z}{H} -\frac{1}{n+3}(\frac{H-z}{H})^{n+3}]$

They say that in order to create Figure 6.2, they use the following values.

  • $H = 100m$
  • $\alpha = 10^{\circ}$
  • $T_s = -10^{\circ}C$
  • $q^{\perp}_{geo} = 50 W m^{-2}$
  • $n = 3$
  • $A = 10^{-16}a^{-1}Pa^{-3} = 10^{-16}a^{-1}kg^{-3} m^{3} s^{6}$
  • $\rho = 910 kg m^{-3}$
  • $\kappa = 2.1 W m^{-1} K^{-1} = 2.1 \mbox{kg } \mbox{ m} \mbox{ s}^{-3}K^{-1}$
  • $g = 9.81 m s^{-2}$

I noticed that when I plug these values into the above equations I don't get the same results as in the figure. I've seen other sources with similar graphs, for example Figure 2 showing $v_x$here, so I know what they did is correct and I feel like I'm missing something.

I checked the units carefully and they line up correctly for the $v_x$ equation. For the $T$ equation I noticed that I need to convert $T_s$ to Kelvin first (ie add 273) before plugging into the equation, and then I needed to convert the final result of that $T$ equation back to Celsius (ie by subtracting 273). This still does not give me a matching graph to what is provided in Figure 6.2.

Why aren't these figures matching when I plug the values into the equations? I've been through the math and coding a bunch and they look ok so I feel like my problem is stemming from a misunderstanding of glacier dynamics. Is there some trick that I'm missing? Thanks!

Here is the original Figure 6.2 enter image description here

Here is what I got enter image description here

Here is the Matlab code I used to get my results in case that helps.

alpha = 10;
A = 10e-16;
rho = 910;
g = 9.81;
kappa = 2.1;
n = 3;
q_geo  = 50;
Ts = -10;

z_grid1 = linspace(0,100,401);
v_x = (2.*A.*(1/(n+1)).*(rho.*g.*sin(alpha)).^n)*(H^(n+1) -(H-z_grid1).^(n+1)); %equation 6.13
title('v_x graph')
xlabel('velocity [m/a]')
ylabel('Height above the base z [m]')

z_grid2 = linspace(0,100,401);
T = (Ts+273) + (q_geo./kappa).*(H-z_grid2)+(2.*A.*H.^(n+3)).* (((rho.*g.*sin(alpha))^(n+1))./ (kappa.*(n+2))).*(1-(z_grid2 ./H) - (1/(n+3)).*((H-z_grid2)./H).^(n+3)); %equation 6.28
title('T graph')
xlabel('T [Degrees Celcius]')
ylabel('Height above the base z [m]')

1 Answer 1


In Matlab, trigonometry functions work with angles in radians, see for example the description of the sin function:

Input angle in radians, specified as a scalar, vector, matrix, or multidimensional array.

So either you keep your code but convert your 10° angle in radians (change to alpha = 0.174533), or you can use the sind function, which also calculates the sine but takes angles in degrees as input.

Then, you also have a factor of 10 error in your definition of $A$: $10^{-16}$ should be noted 1e-16. That fixes the first graph (velocity profile):

enter image description here

The second still looks weird... You also have an error in the definition of $q^{\perp}_{geo}$: the unit is mW, so it should be converted to q_geo=0.05;. This fixes the first line of the equation: if you put the rest in comment, it reproduces well the linear, dashed line, up to a temperature of -7.619°C for $z$ = 0 m. So the issue must be in the second part, although I cannot figure out where.

  • 1
    $\begingroup$ Thank you immensely for your help!! It is very appreciated! $\endgroup$
    – k12345
    Feb 28, 2023 at 13:26
  • 1
    $\begingroup$ You're welcome! I cannot find the issue with the second part of equation 6.28. I am not saying that the source is wrong, but it could be wrong. The full book has an errata correcting several equations (which is perfectly normal, I've never seen a book without any error). Equation 6.28 is not included in the errata, but it could just be that no one has found/reported the error (if any) yet. Keep working on it, and if you feel confident that something is off, you can write to the authors. They are generally glad that someone help them fix their work. $\endgroup$ Feb 28, 2023 at 14:50
  • 1
    $\begingroup$ I appreciate this very much!! Thank you! $\endgroup$
    – k12345
    Feb 28, 2023 at 16:07
  • 1
    $\begingroup$ I have found the issue in 6.28! Just updating this for fun. The A value used in 6.28 has to be in $Pa^{-3} s^{-1}$ in order to have the units match up $\endgroup$
    – k12345
    Mar 3, 2023 at 14:11
  • 1
    $\begingroup$ Makes sense, good catch! You can write this as an answer, it's OK to answer your own question, see earthscience.stackexchange.com/help/self-answer. It would help future readers who encounter the same issue (although I admit that it is unlikely in your case). $\endgroup$ Mar 3, 2023 at 15:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.