Finite difference method not converging to correct steady state or conserving area?

Since a fully NDSolve-based solution is acceptable for you, let me give you one. You simply need the magic of "Pseudospectral" or a dense enough 4th order spatial discretization:

mol[n_, o_:"Pseudospectral"] := {"MethodOfLines", 
  "SpatialDiscretization" -> {"TensorProductGrid", "MaxPoints" -> n, 
    "MinPoints" -> n, "DifferenceOrder" -> o}}

ic[θ_] = Piecewise[{{1.001, -2 Pi/3 < θ < -Pi/3}}, 1.999];

(* Solution 1 *)    
sol = NDSolveValue[{pde, h[-5 π/6, t] == 1.999, h[-π/6, t] == 1.999, 
   h[θ, 0] == ic@θ}, h, {θ, -5 π/6, -π/6}, {t, 0, 10}, Method -> mol[101]]

(* Solution 2, NDSolveValue will spit out eerri warning but it doesn't matter. *)
sol2 = NDSolveValue[{pde, h[-5 π/6, t] == 1.999, h[-π/6, t] == 1.999, 
    h[θ, 0] == ic@θ}, h, {θ, -5 π/6, -π/6}, {t, 0, 10}, 
   Method ->mol[400, 4]];

pic = RegionPlot[1 <= x^2 + y^2 <= 4, {x, -2, 2}, {y, -2, 2}, PlotStyle -> Opacity[0.1]];

ListAnimate@Table[Show[PolarPlot[sol[θ, t], {θ, -5 π/6, -π/6}, 
  PlotRange -> 2, PlotStyle -> Black], pic], {t, 0, 10, 0.1}]

ListAnimate@Table[Show[PolarPlot[sol2[θ, t], {θ, -5 π/6, -π/6}, 
  PlotRange -> 2, PlotStyle -> Black], pic], {t, 0, 10, 0.1}]

OK, now the remaining problem is what's wrong with the FD scheme. Somewhat awkwardly, I found the root of evil is actually a simple mistake: the definition of ds should be

ds = (-(π/6) - -((5 π)/6))/(n - 1)

With this line corrected, your code outputs desired result, too.

Another thing I want to mention is, though (it's surprising to me that) NDSolve still manages to handle implicit vector in your case, it's much faster to define the list of unknown functions explicitly. The following approach is about 180 times faster than yours:

Clear@h;
sol = With[{h = h[#][t] & /@ Range@n}, 
    NDSolveValue[{D[h, t] == dh[h], (h /. t -> 0) == h0}, 
     Head /@ h, {t, 0, 10}]]; // AbsoluteTiming
(* {0.770651, Null} <- on contrast your approach takes about 140 seconds on my laptop*)


valh = Transpose[#["ValuesOnGrid"] & /@ sol];
ti = Flatten@sol[[1]]["Grid"];
s = Array[# &, n, {-5/6 π, -π/6}];

h = ListInterpolation[valh, {ti, s}]; // AbsoluteTiming
(* {0.860209, Null} *)

Introduction

A lack of time to write up an answer ironically provided time to reflect on the problem, and some nagging uncertainties about some issues contributed to the delay.

The slowness of the OP's code can be seen by a simple analysis, which reveals that some expensive calculations are repeated multiple times for each step in the time integration. A more important concern is whether area is conserved by the DE. The OP asserts yes, but I'm thinking not. I wonder which of us is correct. This issue has slowed my response, because I assumed area was conserved, which led to things not making sense. I first thought I could explain the change in area in terms of the finite-difference discretization. It does contribute to it, but not as much as the underlying differential equation. This was confusing to me until I accepted the hypothesis that area is not conserved. Similarly, a difference of opinion about the proper value of the grid spacing ds has made me hesitate some.

Code analysis

The OP's code for dh[] is clearly written and this is plugged into a simple straight-forward NDSolve call, so the whole is not hard to analyze.

Dependence and independence

I have found that I save myself some headaches if variables that depend on one another are constructed in code that reflects the dependency. Let's look at some of the quantities in the code, in particular s and ds.

First, ds and s are intimately linked to each, but their codes are independent. This is a weakness and a potential source for bugs. We should have ds = s[[2]] - s[[1]] (since it is a regular grid). Similarly, h0 depends on n but is not defined in terms of n; indeed, it probably ought to be a function of θ and applied to s.

Control flow

The principal thing to notice is that in dh[], each of d1[#] and d2[#] is repeated several times and each entails a ListCorrelate[] computation to be repeated. This will slow things down unnecessarily. By contrast, @xzczd's fast solution calculates dh[h] just once; even though it is a symbolic computation and somewhat slower than a numerical one, it is not nearly as slow as repeating calculations again and again throughout the steps of the time integration.

Leveraging Mathematica

1. As shown in the tutorial, The Numerical Method of Lines, NDSolve`FiniteDifferenceDerivative[] will compute finite difference derivatives. It is virtually the same as the OP's ListCorrelation[] codes (with the appropriate "DifferenceOrder", with one difference: FiniteDifferenceDerivative[] uses a higher order approximation at the end points of the grid instead of padding (as the OP does with {2.}). The OP's results may be obtained by manually padding in the same way and dropping the first and last results from FiniteDifferenceDerivative[]. For instance, these give the same results as d1 and d2 respectively:

s = Array[# &, n, {-5/6 π, -π/6}];
ds = s[[2]] - s[[1]];  (* == Pi/300 for n = 201 *)
paddedgrid = Join[{First[s] - ds}, s, {Last[s] + ds}];
paddedh0 = Join[{2.}, h0, {2.}];

NDSolve`FiniteDifferenceDerivative[1, paddedgrid,    (* == d1[h0] *)
    "DifferenceOrder" -> 2][paddedh0][[2 ;; -2]]
NDSolve`FiniteDifferenceDerivative[2, paddedgrid,    (* == d2[h0] *)
    "DifferenceOrder" -> 1][paddedh0][[2 ;; -2]]

Note that NDSolve`FiniteDifferenceDerivative does not depend on ds and can be used just in terms of s.

2. A minor point: The following is equivalent to the OP's use of ListCorrelation, padding h with 2. on the left and right:

d1 = ListCorrelate[{-0.5, 0, 0.5}/ds, h, {-2, 2}, 2.]
d2 = ListCorrelate[{1, -2, 1}/ds^2, h, {-2, 2}, 2.]

3. In addition, the OP's use of dh[h[t]], with dh[h_List], means that NDSolve will treat dh[] as a numeric black box, whereas , @xzczd's dh[h] constructs a (large) system of symbolic equations. It seems clear, from the difference in the number of steps shown below, that this affects the integration method used by NDSolve.

Improvement of speed in dh[]

Aside from eliminating the redundant computations outline above, the use of machine numbers, packed arrays, and vectorized operations, as discussed in Performance tuning in Mathematica?, can be used to increase speed. Compile also helps quite a bit. (The code for Compile was constructed from the OP's ode, and I used $ variables to avoid conflicting with other things.) One could do it with the padded grid shown above, which would add a little to the computation time, but I'll do it with the simpler approach.

ClearAll[dh2, s];

n = 201;
s = Developer`ToPackedArray@N@Array[# &, n, {-5/6 π, -π/6}];

dhC = Compile[{{$h, _Real, 1}, {$d1, _Real, 1}, {$d2, _Real, 1},
               {$s, _Real, 1}, {$c, _Real, 1}},
   -2 $d1 (-2 + $h)^2 (-1 + $h) ($c $h + $d1 $s) - 
    2 $d1 (-2 + $h) (-1 + $h)^2 ($c $h + $d1 $s) -
    (-2 + $h)^2 (-1 + $h)^2 (2 $c $d1 + $d2 $s - $h $s)
   ];

With[{ (* compute derivatives and Sin[s], Cos[s] just once *)
   d1FN = NDSolve`FiniteDifferenceDerivative[1, s, "DifferenceOrder" -> 2],
   d2FN = NDSolve`FiniteDifferenceDerivative[2, s, "DifferenceOrder" -> 1],
   sin = Sin[s],
   cos = Cos[s]},
 dh2[h_List] := With[{d1 = d1FN@h, d2 = d2FN@h},
   dhC[h, d1, d2, sin, cos]]
 ];

Performance comparison

I wanted to compare evaluations of the right-hand side of the ODE, steps, and time. I added an option EvaluationMonitor :> ++evals to the NDSolve code for each case. I timed the solutions with code like the following one for my case:

{sol2} = h /. NDSolve[{h'[t] == dh2[h[t]], h[0] == h0}, h, {t, 0, 10}, 
     EvaluationMonitor :> ++evals]; // AbsoluteTiming

Oddly, adding EvaluationMonitor to @xzczd's method increased the time by a huge factor (it took twice as long as mine), so I retimed it without the monitor, where @xzczd's took half the time mine did. These were done with a corrected ds.

One can see that dh2 evaluates about twice as fast, counted per evaluation, than @xzczd's method and 100 times faster than the OP's. Whether that is an entirely fair comparison depends on how much evaluating the ODE accounts for the execution time of NDSolve. It should account for a lot of it, but that's just my guess.

Area preserving?

I figure I should share this. If I made a mistake, well, someone can point it out. If not, then, maybe it will help the OP. If we consider what it means to be area-preserving, we should consider the integral $$A = \int_{-\pi}^\pi \text{$1\over2$}\, h(t,\theta)^2\;d\theta \,,$$ where $h$ is extended to $-\pi \le \theta \le \pi$ by defining $h(t,\theta) = 2$ (that is, the greater radius) outside $-5\pi/6 \le \theta \le -\pi/6$. Then, by my lights, differentiating under the integral sign, the rate of change of the area is given by $${dA \over dt} = \int_{-\pi}^\pi h(t,\theta)h_t(t,\theta)\;d\theta$$ and $h_t$ can be written in terms of $h$, $h_\theta$, $h_{\theta\theta}$ using the PDE.

Of course, an equilibrium solution will preserve area. But if we perturb an equilibrium solution, then $dA/dt$ is nonzero:

ClearAll[dhdt, dAdt];
Block[{h, θ},
 dhdt[h_, θ_] = -D[(h[θ] - 2)^2 (h[θ] - 1)^2 D[h[θ] Sin[θ], θ], θ];
 ];
dAdt[hfn_] := NIntegrate[hfn[θ] dhdt[hfn, θ], {θ, -Pi, Pi}];

(* background *)
θ1 = -2.4539664700963293`;
background = Show[
   RegionPlot[1 <= x^2 + y^2 <= 4, {x, -2, 2}, {y, -2, 2}, 
    PlotStyle -> Opacity[0.1], AspectRatio -> 1],
   PolarPlot[
    2 Sin[θ1]/Sin[θ], {θ, θ1, -(Pi + θ1)}, 
    PlotStyle -> Black]
   ];

solPlot[hfn_] := Show[
   background,
   PolarPlot[hfn, {θ, -5/6 π, -π/6},
    ColorFunction -> (ColorData["RedGreenSplit"][
        0.5 (1 + Clip[100 dhL /. h -> Function @@ {θ, hfn} /. θ -> #3, {-1,1}])] &),
    ColorFunctionScaling -> False,
    PlotRange -> {{-2, 2}, {-2, 2}}]
   ];

Block[{Print, θ2 = θ1, eps = +0.1, hf},
 hf = Piecewise[{
     {2 Sin[θ2] Csc[#] + eps (# - θ2) (# + Pi + θ2), θ2 < # < -Pi - θ2}
     }, 2] &;
 Show[solPlot[hf[θ]], PlotLabel -> Row[{"dA/dt = ", dAdt[hf]}]]
 ]

(Red means the radius h is decreasing; green means it is increasing (with respect to time in the DE.)

The effect of finite differences and the error in ds

Using Block and the OP's definition of dh, we can see what happens when we use different values for ds. There is some controversy about the correct value for ds, and it is my view that it should be equal to the differences in s. The padding by 2 on either side of h in for ListCorrelate[] in the OP's dh[] corresponds to using s == -5 Pi/ 6 - ds and s == -Pi / 6 + ds and corresponds to the (unused) condition that h is identically equal to 2 outside the interval {-5 Pi/6, -Pi/6}.

If we plot dh[h[s]], where h is the equilibrium solution of the OP, h[s] = 2 Sin[θ1] Csc[s] for a given θ1, and s is the input grid as above, we will see how close to zero the value of dh is for the finite difference scheme.

ClearAll[dhEquilPlot];
SetAttributes[dhEquilPlot, Listable];
θ1 = -2.4539664700963293`;
dhEquilPlot[ds0_] := Block[{ds = ds0},
  ListPlot[
   dh[Clip[2 Sin[θ1] Csc[s], {0., 2.}]],
   PlotLabel -> Row[{"ds = ", ds}], Frame -> True, 
   FrameLabel -> {HoldForm[s], HoldForm[h'[t]]}, DataRange -> MinMax[s], 
   BaseStyle -> {FractionBoxOptions -> {Beveled -> True}}
   ]
  ]

dhEquilPlot /@ {
   {1/n},                               (* OP     (UL) *)
   {(-(π/6) - -((5 π)/6))/(n + 2 - 1),  (* xzczd  (LL) *)
    (-(π/6) - -((5 π)/6))/(n - 1)}      (* M E2   (LR) *)
   } // GraphicsGrid

We can see how much closer to zero dh is for ds = Pi/300, which is also the grid spacing in s. It is also clear that in all three cases, the radius h is pulled toward the center (dh < 0) near where the level of the fluid contacts the outside circle (where h reaches 2). Of course the actual evolution from the OP's starting configuration will be slightly different. But since physically, the model should be area-preserving, there is little point in exploring this beyond determining the value of ds until the issue with area is cleared up.