So I am trying to solve the movement in space and time of a spreading gravity current. The interface satisfies the following PDE:

$ \frac{\partial h}{\partial t} = \frac{\partial}{\partial x}\left(h^3 \frac{\partial h}{\partial x}\right) $.

The nose, $h(x,t) = 0$, of the current can move forward in space. If we say the initial condition is $h(x,0) = 1-x$, then the nose, $x_N$ is initially at $x = 1$. The base of the current can be fixed at $h(0,t) = 1$. So I have one IC and one BC, I need another BC to close the system. Now, as time moves forward the nose (fixed at h = 0) propagates according to the expression:

$\int_0^{x_N} h \, dx = \frac{1}{2} + t$

or kinematically $\dot{x_N} = -h^2 h_x$.

Clearly the second boundary condition in $x$ needs to come from the nose condition, but I have no idea how to input an integral boundary condition or the kinematic condition. Secondly, I have no idea how to deal with xmax in NDSolve as clearly my xmax is moving forward with each time step.

Here is my attempt, which gives a "There are fewer dependent variables error"

NDSolve[{D[h[x, t], t] == D[h[x, t]^3 D[h[x, t], x], x], h[0, t] == 1,
h[x, 0] == 1 - x (D[h[x, t], t] /. x -> 1) == -h[x, t]^2 (D[h[x, t], x] /. 
  x -> 1)}, h, {x, 0, 1}, {t, 0, 1}]


So I solved this problem by amending the code in this Stefan problem post. I modified the statement of the problem slightly to work more closely to @ybeltukov excellent solution, and actually solved the following system:

$ \frac{\partial h}{\partial t} = \frac{\partial}{\partial x}\left(h^3 \frac{\partial h}{\partial x}\right) \\ \dot{s} = -h^2 h_x \\ s(0) = 0 \\ h(x,0) = \begin{cases}1 \quad \text{for } x = 0 \\ 0 \quad \text{otherwise}\end{cases} \\ h(s(t),t) = 0 \\ h(0,t) = 1 $.

I made the same change to a normalised variable for my PDE and then amended finite difference code:

n = 100;
\[Delta]\[Xi] = 1./n;

ClearAll[dv, t];
dv[v_List] := 
With[{s = First@v, u = Rest@v}, 
With[{ds = u[[-1]]^3/(s \[Delta]\[Xi]), \[Xi] = N@Range[n - 1]/n, 
 d1 = ListCorrelate[{-0.5, 0, 0.5}/\[Delta]\[Xi], #] &, 
 d2 = ListCorrelate[{1, -2, 1}/\[Delta]\[Xi]^2, #] &}, 
Prepend[3 u^2 d1[#]^2/s^2 + u^3 d2[#]/s^2 + \[Xi] ds d1[#]/s &@
  Join[{1}, u, {0.}], ds]]];
 s0 = 0.001;
 v0 = Flatten@Prepend[ConstantArray[0.001, n - 2], {s0, 1.}];
 sol = NDSolve[{v'[t] == dv[v[t]], v[0] == v0}, v, {t, 0, 50}][[1, 1, 

Returning from the normalised variable is identical code.

You can then make a nice time series with this code:

Table[Plot[u[t, x], {x, 0, 10}, PlotRange -> {{0, 10}, {0, 3}}], {t, 
 0, 50, 0.1}]]

Which returns a spreading current as desired.

Spreading viscous gravity current

  • $\begingroup$ Could you explain what do you mean by the kinematic condition $\dot{x}_N=-h^2h_x$? Do I understand correctly that the right-hand-side should be evaluated at $x=0$? or at $x=x_N$ (for which it vanishes)? $\endgroup$
    – yohbs
    Commented May 12, 2016 at 19:30
  • $\begingroup$ Also, your initial condition is inconsistent with the integral constraint because at $t=0$ you have $$\int_0^{x_N} h\, dx=\int_0^1(1-x)\,dx=\frac{1}{2}\ne t=0$$ $\endgroup$
    – yohbs
    Commented May 12, 2016 at 19:38
  • $\begingroup$ Similar to mathematica.stackexchange.com/questions/109201/… $\endgroup$ Commented May 12, 2016 at 20:28
  • $\begingroup$ @ojlm I'm not sure if your statements about the moving nose are correct. You might wish to consult my solution, and particularly the steady state solution. $\endgroup$ Commented May 12, 2016 at 20:35
  • $\begingroup$ @Dr.WolfgangHintze The nose needs to propagate downstream as this is a model for a fluid mechanical process. Pinning the interface to the bottom boundary at x=1 is not physically what happens. I have given a link to a paper in my other comment. Governing equations are on page 3. I have simplified the propagation flux condition as I am only considering a one layered spreading current in this model. $\endgroup$
    – mch56
    Commented May 12, 2016 at 20:44

1 Answer 1


The PDE can be easily solved numerically if a simple additional boundary condition is imposed.

Adding a second bundary condition h[1,t] = 0 we find numerically, without any error message

hh[x_, t_] = 
  h[x, t] /. 
   NDSolve[{D[h[x, t], t] == D[h[x, t]^3 D[h[x, t], x], x], h[0, t] == 1, 
      h[1, t] == 0, h[x, 0] == (1 - x)}, h[x, t], {x, 0, 1}, {t, 0, 5}][[1]];

Plot3D[hh[x, t], {x, 0, 1}, {t, 0, 3}, PlotRange -> {0, 1}, 
 PlotLabel -> "Solution of a PDE", AxesLabel -> {"x", "t", "h[x,t]"}]

enter image description here

Executing NDSolve without a secod boundary condition leads to about the same graph and the error message that a second boundary condition is missing.

The time dependence at various positions is given by the graph

Plot[{hh[0.2, t], hh[0.5, t], hh[0.8, t], hh[0.95, t]}, {t, 0, 3}, 
 PlotRange -> {0, 1}, 
 PlotLabel -> 
  "Solution of a PDE\nTime dependence for various locations x = 1 - h[x,0]", 
 AxesLabel -> {"t", "h[x,t]"}]

enter image description here

The asymptotic function for large t is given by

ha[x_] = (1 - x)^(1/4)

which solves the steady state ODL derived from the original PDE.

We could study the approch of the solution to the equlibirum by linearizing the equation, i.e. letting

h = hL[x] + ha[x]

and dropping all non linear terms in hL.

Partial solutions are given by Bessel functions. The interested reader will surely be able to follow this path on his own.

Moving structures

Let us have a look at moving structures allowed by the PDE.


h[x,t] = g[x - v t]

The PDE reduces to an ODE which can be solved

gg[x_] = g[x] /. 
   DSolve[g'[x] v == 3 g[x]^2 g'[x]^2 + g[x]^3 g''[x], g[x], x][[1]] /. 
  C[2] -> 0

(* Out[341]= InverseFunction[-((C[1]^3 Log[C[1] + v #1])/v^4) + (C[1]^2 #1)/v^3 - (
    C[1] #1^2)/(2 v^2) + #1^3/(3 v) &][x] *)


h1[s_] = (gg[x][[0, 1]] /. {C[1] -> -10, v -> 1} // FullSimplify)@s

(* Out[393]= 100 s + 5 s^2 + s^3/3 + 1000 Log[-10 + s] *)


Plot[Re[h1[s]], {s, 0, 10}, PlotRange -> {0, 3540}, 
 PlotLabel -> "'Nose' profile", AxesLabel -> {"x - t", "g"}]

enter image description here

This section is just a sketch which must be worked out more thoroughly.

  • $\begingroup$ Thanks for the response. I know it's possible to solve the PDE with h[1,t] = 0 boundary condition. However, since this is the propagation of a current, the nose needs to move over the bottom boundary. I want to be able to to create a time series such that the current moves along down the x axis. If you see this paper, you can see how the current propgates down the x-axis as time progresses. Having the h[1,t] = 0 boundary condition does not allow me this. repository.cam.ac.uk/bitstream/handle/1810/247602/… $\endgroup$
    – mch56
    Commented May 12, 2016 at 20:41
  • $\begingroup$ I appreciate your further analysis about travelling wave solutions. This problem is well known and the best way to study it analytically is via a similarity solution. However, the point is I am using this model in an illustrative way. I am studying a much harder system, that has an interface moving transiently along a boundary. The key thing I am interested in do is to create a time series of the spreading thinning profile as seen in the paper I sent. Would you be able to help me in achieving this? I want to be have the profile spreading and thinning with accordance to the volume BC $\endgroup$
    – mch56
    Commented May 13, 2016 at 8:17

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.