Since you're working over the domain Reals, use Sqrt[#.#] & or better in this case, its square, #.# &.

{pM = {xM, yM, 2}, pD = {0, 0, 0}, pN = {0, 2, h}, pA = {2, 0, 0}, 
  relcir = (xM - 1)^2 + (yM - 1)^2 == 1, pQ = {a, b, c}};
slq = Solve[{pQ - pA == \[Lambda] (pN - pA)}][[1]];
pQ = pQ /. slq;
Block[{Norm = # . # &},
 AbsoluteTiming@
  Reduce[{Norm[pQ - pA] == Norm[pQ - pN] == Norm[pQ - pM] == 
     Norm[pQ - pD], relcir, 0 <= h <= 2, 
    0 <= \[Lambda] <= 1}, {\[Lambda], h}, {xM, yM}, Reals]
 ]
(*  {0.051348, \[Lambda] == 1/2 && h == 3/2}  *)

I suppose I should have mentioned that what guided me was old-Mathematica, when polynomial systems and systems that could be converted to polynomial systems were the principal strength of Mathematica in symbolic equation-solving. The complex (numbers) function Abs[] that appears in Norm[] adds a level of complexity that apparently invokes a lot of algebra or a lot of checking, or both. While equation-solving has greatly improved over the years, this drawback of using Norm[] in real-algebraic equations has persisted. Similarly, Sqrt[] means extra branch-checking and extra time. The Reals specification means that Reduce[] has to check that the value of each Sqrt[] is real in its operations.