<p>In Matlab, trigonometry functions work with angles in radians, see for example the <a href="https://www.mathworks.com/help/matlab/ref/sin.html" rel="nofollow noreferrer">description</a> of the <code>sin</code> function:</p>
<blockquote>
<p>Input angle in radians, specified as a scalar, vector, matrix, or
multidimensional array.</p>
</blockquote>
<p>So either you keep your code but convert your 10° angle in radians (change to <code>alpha = 0.174533</code>), or you can use the <code>sind</code> function, which also calculates the sine but takes angles <a href="https://www.mathworks.com/help/matlab/ref/sind.html#btg_ljb-3" rel="nofollow noreferrer">in degrees</a> as input.</p>
<p>Then, you also have a factor of 10 error in your definition of <span class="math-container">$A$</span>: <span class="math-container">$10^{-16}$</span> should be noted <code>1e-16</code>. That fixes the first graph (velocity profile):</p>
<p><a href="https://i.sstatic.net/pJdLH.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/pJdLH.png" alt="enter image description here" /></a></p>
<p>The second still looks weird... You also have an error in the definition of <span class="math-container">$q^{\perp}_{geo}$</span>: the unit is mW, so it should be converted to <code>q_geo=0.05;</code>. 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 <span class="math-container">$z$</span> = 0 m. So the issue must be in the second part, although I cannot figure out where.</p>

CC BY-SA 4.0
2023-02-28T13:26:08.200Z
Thank you immensely for your help!! It is very appreciated!
k12345
CC BY-SA 4.0
2023-02-28T14:50:04.360Z
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](https://zenodo.org/record/3993753) 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.
Jean-Marie Prival
CC BY-SA 4.0
2023-02-28T16:07:34.053Z
I appreciate this very much!! Thank you!
k12345
CC BY-SA 4.0
2023-03-03T14:11:06.840Z
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
k12345
CC BY-SA 4.0
2023-03-03T15:55:27.317Z
Makes sense, good catch! You can write this as an answer, it's OK to answer your own question, see https://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).
Jean-Marie Prival
CC BY-SA 4.0
2023-02-28T09:12:27.420Z
2023-02-28T15:57:37.677Z
2023-02-28T15:57:37.677Z
Jean-Marie Prival
<p>I'm trying to recreate the graphs in Figure 6.2 from this <a href="https://ocw.hokudai.ac.jp/wp-content/uploads/2016/03/DynamicsOfIce-2010-Note-all.pdf" rel="nofollow noreferrer">source</a> (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).</p>
<p><span class="math-container">$v_x = \frac{2 A (\rho g sin(\alpha))^{n}}{n+1}(H^{n+1} -(H-z)^{n+1})$</span><br />
<span class="math-container">$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}]$</span></p>
<p>They say that in order to create Figure 6.2, they use the following values.</p>
<ul>
<li><span class="math-container">$H = 100m$</span></li>
<li><span class="math-container">$\alpha = 10^{\circ}$</span></li>
<li><span class="math-container">$T_s = -10^{\circ}C$</span></li>
<li><span class="math-container">$q^{\perp}_{geo} = 50 W m^{-2}$</span></li>
<li><span class="math-container">$n = 3$</span></li>
<li><span class="math-container">$A = 10^{-16}a^{-1}Pa^{-3} = 10^{-16}a^{-1}kg^{-3} m^{3} s^{6}$</span></li>
<li><span class="math-container">$\rho = 910 kg m^{-3}$</span></li>
<li><span class="math-container">$\kappa = 2.1 W m^{-1} K^{-1} = 2.1 \mbox{kg }  \mbox{ m} \mbox{ s}^{-3}K^{-1}$</span></li>
<li><span class="math-container">$g = 9.81 m s^{-2}$</span></li>
</ul>
<p>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 <span class="math-container">$v_x$</span><a href="https://glaciers.gi.alaska.edu/sites/default/files/Notes_icesheetmod_Bueler2014.pdf" rel="nofollow noreferrer">here</a>, so I know what they did is correct and I feel like I'm missing something.</p>
<p>I checked the units carefully and they line up correctly for the <span class="math-container">$v_x$</span> equation. For the <span class="math-container">$T$</span> equation I noticed that I need to convert <span class="math-container">$T_s$</span> to Kelvin first (ie add 273) before plugging into the equation, and then I needed to convert the final result of that <span class="math-container">$T$</span> 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.</p>
<p><strong>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!</strong></p>
<p>Here is the original Figure 6.2
<a href="https://i.sstatic.net/NNHTp.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/NNHTp.png" alt="enter image description here" /></a></p>
<p>Here is what I got
<a href="https://i.sstatic.net/HO8jG.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/HO8jG.png" alt="enter image description here" /></a></p>
<p>Here is the Matlab code I used to get my results in case that helps.</p>
<pre><code>H=100;
alpha = 10;
A = 10e-16;
rho = 910;
g = 9.81;
kappa = 2.1;
n = 3;
q_geo  = 50;
Ts = -10;

figure()
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
plot(v_x,z_grid1)
title('v_x graph')
xlabel('velocity [m/a]')
ylabel('Height above the base z [m]')



figure()
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
plot(T-273,z_grid2)
title('T graph')
xlabel('T [Degrees Celcius]')
ylabel('Height above the base z [m]')
</code></pre>

<glaciology><numerical-modelling><matlab>
Understanding Graphs for Velocity and Temperature Solutions for Parallel Sided Slab on Slope Problem
k12345
# Understanding Graphs for Velocity and Temperature Solutions for Parallel Sided Slab on Slope Problem

(Asked by: username_0 on 2023-02-27)

I'm trying to recreate the graphs in Figure 6.2 from this [source](https://ocw.hokudai.ac.jp/wp-content/uploads/2016/03/DynamicsOfIce-2010-Note-all.pdf) (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](https://glaciers.gi.alaska.edu/sites/default/files/Notes_icesheetmod_Bueler2014.pdf), 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](https://i.sstatic.net/NNHTp.png)](https://i.sstatic.net/NNHTp.png)

Here is what I got
[![enter image description here](https://i.sstatic.net/HO8jG.png)](https://i.sstatic.net/HO8jG.png)

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

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

figure()
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
plot(v_x,z_grid1)
title('v_x graph')
xlabel('velocity [m/a]')
ylabel('Height above the base z [m]')

figure()
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
plot(T-273,z_grid2)
title('T graph')
xlabel('T [Degrees Celcius]')
ylabel('Height above the base z [m]')

```

**Tags:** glaciology, numerical-modelling, matlab


---

## 1 answers

### Answer by username_1 on 2023-02-28. Score: 2 | Accepted answer

In Matlab, trigonometry functions work with angles in radians, see for example the [description](https://www.mathworks.com/help/matlab/ref/sin.html) 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](https://www.mathworks.com/help/matlab/ref/sind.html#btg_ljb-3) 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](https://i.sstatic.net/pJdLH.png)](https://i.sstatic.net/pJdLH.png)

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.

> **Comments:**
> **username_0**: Thank you immensely for your help!! It is very appreciated!  
> **username_1**: 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](https://zenodo.org/record/3993753) 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.  
> **username_0**: I appreciate this very much!! Thank you!  
> **username_0**: 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  
> **username_1**: Makes sense, good catch! You can write this as an answer, it's OK to answer your own question, see https://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).  

---

