NDSolve not handling discontinuity in initial condition

Question

I am trying to solve a PDE in two variables z and t. This PDE is from "Universality of Capillary Rise" by Zhou and Doi (2020). There is additonally a parameter n that can take positive integer values. The initial condition is defined as a piecewise function with a linear component and a zero component. However, when I attempt to use NDSolve on this problem and examine the output at t=0 the output is not close to the linear profile I specified. Further, while my code performs ok for n=1 it does not come close to evaluating correctly for n=2. Rather, for that case it appears that the code gets stuck and never advances the solution in time. Any help would be appreciated!

geqn[n_] := 
  Simplify[
   1/((3 n + 1 ) G[z, t]^n)* 
    D[G[z, t]^(3 n + 1) (1 + 1/G[z, t]^(n + 1) D[G[z, t], z]), z]];
equilBottom[n_] := (.1*n)^(-1/n);
initialEqn[z_, Z0_, n_] := 
  Piecewise[{{equilBottom[n] (1 - (z - .1)/(Z0 - .1)), 
     0.1 <= z <= Z0}, {0, Z0 < z}}];
eqnList[n_] := {D[G[z, t], t] == geqn[n], G[0.1, t] == equilBottom[n],
   G[20, t] == 0, G[z, 0] == initialEqn[z, 1, n]}
sol1 = NDSolve[eqnList[1], G[z, t], {z, .1, 20}, {t, 0, 1000}] // 
   First;
sol2 = NDSolve[eqnList[2], G[z, t], {z, .1, 20}, {t, 0, 50000}] // 
   First;
sol11 = sol1[[1, 2]];
sol21 = sol2[[1, 2]];

The output of these two solutions at various values of t is shown below. The output is along the x-axis and the y-axis is z. Note that the black line is the equilibrium solution of the PDE.

Further, the initial condition in the case of n=1 is shown below. Note the similarity to the output of the second plot while not obeying the initial condition as specificed in the equaiton.

Here is the code for the graphics:

Import["https://pastebin.com/raw/uUJ4rmc1", "Package"]
colors = Table[parulaNew[i], {i, 0, 1, 1/(5 - 1)}];
plots = MapThread[
   ParametricPlot[{sol11 /. {z -> z0, t -> #1}, z0}, {z0, .1, 20}, 
     PlotRange -> {{0, 1}, {0, 20}}, AspectRatio -> 9/16, 
     PlotStyle -> #2] &, {{200, 400, 600, 800, 1000}, colors}];
equiPlot = ParametricPlot[{1/z, z}, {z, .1, 20}, PlotStyle -> Black];
plots = Join[plots, {equiPlot}];
p1 = Show[plots]
plots2 = MapThread[
   ParametricPlot[{sol21 /. {z -> z0, t -> #1}, z0}, {z0, .1, 20}, 
     PlotRange -> {{0, 1}, {0, 20}}, AspectRatio -> 9/16, 
     PlotStyle -> #2] &, {{1000, 2000, 3000, 4000, 50000}, colors}];
equiPlot2 = 
  ParametricPlot[{(2 z)^(-1/2), z}, {z, .1, 20}, PlotStyle -> Black];
plots2 = Join[plots2, {equiPlot2}];
p2 = Show[plots2]
p3 = ParametricPlot[{sol11 /. {z -> z0, t -> 0}, z0}, {z0, .1, 20}, 
  PlotRange -> {{-.1, 1}, {0, 5}}, AspectRatio -> 9/16]