We can solve this problem using the Euler wavelets collocation method and method of lines. Supposed that h[x,t] is a real function, we express Abs[hx] as Sqrt[hx^2], so we have
ClearAll[h, x, n, t]
hx = D[h[x, t], x];
ht = D[h[x, t], t];
n = 1;
sqrtTerm = Sqrt[1 + 4nSqrt[hx^2]];
pde = Sign[-hx]1/2D[h[x, t]*(sqrtTerm - 1), x] == -ht;
Nevertheless it can't be solved with NDSolve due to some reason. Therefore we use wavelets as a method to discretize the equation on x
OEm[m_, x_] :=
Sqrt[2 m +
1] Sum[(-1)^(m - k) x^k Binomial[m, k] Binomial[m + k, k], {k, 0,
m}]; UE[m_, t_] := OEm[m, t];
psi[k_, n_, m_, t_] :=
Piecewise[{{2^((k - 1)/2) UE[m, 2^(k - 1) t - n + 1], (n - 1)/
2^(k - 1) <= t < n/2^(k - 1)}, {0, True}}];
PsiE[k_, M_, t_] :=
Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]];
k0 = 2; M0 = 4; With[{k = k0, M = M0},
nn = Length[Flatten[Table[1, {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]]];
dx = 1/(nn); xl = Table[l*dx, {l, 0, nn}]; xcol =
Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, nn + 1}]; Psijk =
With[{k = k0, M = M0}, PsiE[k, M, t1]]; Int1 =
With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Int2 = Integrate[Int1, t1];
Psi[y_] := Psijk /. t1 -> y; int1[y_] := Int1 /. t1 -> y;
int2[y_] := Int2 /. t1 -> y;
wA = Table[wa[i][t], {i, nn}]; wB = Table[wb[i][t], {i, 2}];
w2[x_] := wA . Psi[x]; w1[x_] := wA . int1[x] + wB[[1]];
w0[x_] := wA . int2[x] + wB[[1]] x + wB[[2]];
eqw = With[{w =
w0[x]}, -D[w, t] == (-(1/2))*
Sign[Derivative[1, 0][h][x,
t]]*(Derivative[1, 0][h][x,
t]*(-1 +
Sqrt[1 + 4*Sqrt[Derivative[1, 0][h][x, t]^2]]) + (2*
h[x, t]*Sign[Derivative[1, 0][h][x, t]]*
Derivative[2, 0][h][x, t])/
Sqrt[1 + 4*Sqrt[Derivative[1, 0][h][x, t]^2]]) /. {h[x, t] ->
w0[x], Derivative[1, 0][h][x, t] -> w1[x],
Derivative[2, 0][h][x, t] -> w2[x]}];
eqnw = Table[eqw, {x, xcol}];
icx = With[{w = w0[x]}, w == 1 /. t -> 10^-3];
ic = Table[icx, {x, xcol}];
bc = With[{w = w0[x]},
Join[{w1[x] == 0} /. x -> 0, {w == 0} /. x -> 1]]; varAll =
Join[wA, wB];
icn = Join[ic, bc /. t -> 10^-3]; eqn = Join[eqnw, D[bc, t]]; var1 =
D[varAll, t];
{vec, mat} = CoefficientArrays[eqn, var1];
f = Inverse[mat // N] . (-vec); vr0 = varAll /. t -> 10^-3;
{v0, mat0} = CoefficientArrays[icn, vr0];
sol0 = LinearSolve[mat0, -v0];
icn0 = Table[vr0[[i]] == sol0[[i]], {i, Length[vr0]}];
sol1 = NDSolve[{Table[var1[[i]] == f[[i]], {i, Length[var1]}], icn0},
varAll, {t, 10^-3, 10}];
Visualization
ParametricPlot3D[{t, x, w0[x] /. sol1[[1]]}, {t, 10^-3, 10}, {x, 0,
1}, ColorFunction -> Hue, AxesLabel -> {"t", "x", "h"},
Boxed -> False, BoxRatios -> {1, 1, 1}, PlotPoints -> {40, 35},
MaxRecursion -> 1, MeshStyle -> {Red, Blue}]
Update 1. This problem can be solved with FDM as follows
L = 1; dx = L/20; xgrid = Range[0, L, dx]; nn = Length[xgrid];
M2 = NDSolve`FiniteDifferenceDerivative[Derivative[2], xgrid, DifferenceOrder -> 2]["DifferentiationMatrix"];
M1 = NDSolve`FiniteDifferenceDerivative[Derivative[1], xgrid, DifferenceOrder -> 2]["DifferentiationMatrix"];
wA = Table[wa[i][t], {i, nn}];
w1 = M1 . wA; w2 = M2 . wA;
rhs = (-(1/2))*Sign[Derivative[1, 0][h][x, t]]*Derivative[1, 0][h][x, t]*(-1 + Sqrt[1 + 4*Sqrt[Derivative[1, 0][h][x, t]^2]]) -
(h[x, t]*Derivative[2, 0][h][x, t])/Sqrt[1 + 4*Sqrt[Derivative[1, 0][h][x, t]^2]] /. {h[x, t] -> wA, Derivative[1, 0][h][x, t] -> M1 . wA,
Derivative[2, 0][h][x, t] -> M2 . wA};
eq = Table[D[wA[[i]], t] == -rhs[[i]], {i, 2, nn - 1}];
bc = {(M1 . wA)[[1]] == 0, wA[[-1]] == 0};
ic = Table[wA[[i]] == 1 /. t -> 10^(-3), {i, 2, nn - 1}]; icn = Join[ic, bc /. t -> 10^(-3)];
eqn = Join[eq, D[bc, t]]; var1 = D[wA, t];
{vec, mat} = CoefficientArrays[eqn, var1];
ff = Inverse[N[mat]] . (-vec); vr0 = wA /. t -> 10^(-3);
{v0, mat0} = CoefficientArrays[icn, vr0];
sol0 = LinearSolve[mat0, -v0];
icn0 = Table[vr0[[i]] == sol0[[i]], {i, Length[vr0]}];
Off[General::partd];
f[t_, x_] := Evaluate[ff /. Table[wa[i][t] -> x[[i]], {i, nn}]]
rk2[f_, h_][{t_, x_}] := Module[{k1, k2}, k1 = f[t, x];
k2 = f[t + h/2, x + h k1/2];
{t + h, x + h k2}];
tf = 2; dt = 1/1000; sol =
NestList[rk2[f, dt], {0, icn0[[All, 2]]}, Round[tf/dt]];
Visualization
h = Interpolation[
Flatten[Table[{{sol[[j, 1]], xgrid[[i]]}, sol[[j, 2, i]]}, {j,
Length[sol]}, {i, nn}], 1], InterpolationOrder -> 1]
Plot3D[h[t, x], {t, 0, 2}, {x, 0, 1}, ColorFunction -> Hue,
PlotTheme -> "Scientific", AxesLabel -> {"t", "x", "h"},
PlotRange -> All, PlotPoints -> 50, MeshStyle -> {Red, Blue}]
Note, in this example we used the Runge-Kutta second order algorithm on $0\le t \le 2$ with step $10^{-3}$ and difference scheme of second order on $0\le x\le 1$ with step 1/20.
Update 2
This is more advanced code with the Runge-Kutta 4th order algorithm on $0\le t\le 10$ and difference scheme of 4th order on $0\le x\le1$
L = 1; tend = 10; dx = L/20; xgrid = Range[0, L, dx]; nn = Length[xgrid];
M2 = NDSolve`FiniteDifferenceDerivative[Derivative[2], xgrid,
DifferenceOrder -> 4]@"DifferentiationMatrix"; M1 =
NDSolve`FiniteDifferenceDerivative[Derivative[1], xgrid,
DifferenceOrder -> 4]@"DifferentiationMatrix";
wA = Table[wa[i][t], {i, nn}];
w1 = M1 . wA; w2 = M2 . wA;
rhs = (-(1/2))Sign[Derivative[1,0][h][x,t]](Derivative[1,0][h][x,t](-1+Sqrt[1+4Sqrt[Derivative[1,0][h][x,t]^2]])+(2h[x,t]Sign[Derivative[1,0][h][x,t]]Derivative[2,0][h][x,t])/Sqrt[1+4Sqrt[Derivative[1,0][h][x,t]^2]]) /. {h[x, t] -> wA, Derivative[1, 0][h][x, t] -> M1 . wA,
Derivative[2, 0][h][x, t] -> M2 . wA};
eq = Table[D[wA[[i]], t] == -rhs[[i]], {i, 2, nn - 1}];
bc = {(M1 . wA)[[1]] == 0, wA[[-1]] == 0};
ic = Table[wA[[i]] == 1 /. t -> 10^-3, {i, 2, nn - 1}]; icn =
Join[ic, bc /. t -> 10^-3];
eqn = Join[eq, D[bc, t]]; var1 = D[wA, t];
{vec, mat} = CoefficientArrays[eqn, var1];
ff = Inverse[mat // N] . (-vec); vr0 = wA /. t -> 10^-3;
{v0, mat0} = CoefficientArrays[icn, vr0];
sol0 = LinearSolve[mat0, -v0];
icn0 = Table[vr0[[i]] == sol0[[i]], {i, Length[vr0]}];
Off[General::partd];
f[t_, x_] := Evaluate[ff /. Table[wa[i][t] -> x[[i]], {i, nn}]]
rk4[f_, h_][{t_, y_}] := Module[{k1, k2, k3, k4}, k1 = f[t, y];
k2 = f[t + h/2, y + h k1/2];
k3 = f[t + h/2, y + h k2/2];
k4 = f[t + h, y + h k3];
{t + h, y + h/6*(k1 + 2 k2 + 2 k3 + k4)}]
Example of usage
dt1 = 1/1000; tf = 10; sol1 =
NestList[rk4[f, dt1], {0, icn0[[All, 2]]}, Round[tf/dt1]];
Visualization
h = Interpolation[
Flatten[Table[{{sol1[[j, 1]], xgrid[[i]]}, sol1[[j, 2, i]]}, {j,
Length[sol1]}, {i, nn}], 1], InterpolationOrder -> 8];
Plot[Table[h[t, x], {t, {.1, .5, 1, 2, 5, 10}}] // Evaluate, {x, 0,
1}, Frame -> True, FrameLabel -> {"x", "H"},
PlotLegends -> Table[Row[{"t =", t}], {t, {.1, .5, 1, 2, 5, 10}}]]
