Python3, 523 bytes

E=enumerate
M=[(1,0),(-1,0),(0,-1),(0,1),(1,1),(-1,1),(-1,-1),(1,-1)]
def f(L):
 d,t={(x,y):v for x,r in E(L)for y,v in E(r)},[]
 for i in d.values():
  q={j for j in d if d[j]==i and j not in t}
  while q:
   v=[*q][0];Q,s={v},{v}
   while Q:e=[*Q][0];Q-={e};x,y=e;l={V for X,Y in M if(V:=(x+X,y+Y))in d and d[V]<=d[v]and V not in s};Q={*Q,*l};s={*s,*l}
   q-=s;t+=[*s]*all(d.get((x+X,y+Y),-1)>d[(x,y)]or(x+X,y+Y)in[*t,*s]for x,y in s for X,Y in M)
 for x,y in t:L[x][y]='*'
 return'\n'.join(''.join(map(str,i))for i in L)

Try it online!

Python, 246 chars

import os
a=list(os.read(0,2e3))
w=a.index('\n')+1
a+=' '*w
def f(p,t):
    if e<a[p]or p in t:return
    t[p]=1
    return'*'>a[p]or any(f(p+d,t)for d in(~w,-w,-w+1,-1,1,w-1,w,w+1))
z=0
for e in a:
    if(' '<e)*~-f(z,{}):a[z]='*'
    z+=1
print''.join(a[:~w])

The solution works by doing a DFS from each position to determine whether or not to fill.

If trailing whitespace on each line is allowed, it may be shortened by using w=80 and padding the input lines with whitespace to 80 chars.

Python, 478 471 characters

(Not including comments. 452 450 characters not including the imports.)

import sys,itertools
i=[list(x)for x in sys.stdin.read().strip().split('\n')]
h=len(i)
w=len(i[0])
n=h*w
b=n+1
e=range(h)
d=range(w)
# j is, at first, the adjancency matrix of the graph.
# The last vertex in j is the "drain" vertex.
j=[[[b,1][(t-r)**2+(v-c)**2<=1 and i[r][c]>=i[t][v]] for t in e for v in d]+[[b,1][max([r==0,r>h-2,c==0,c>w-2])]]for r in e for c in d]+[[0]*b]
r=range(b)
for k,l,m in itertools.product(r,repeat=3):
    # This is the Floyd-Warshall algorithm
    if j[l][k]+j[k][m]<j[l][m]:
        j[l][m]=j[l][k]+j[k][m]
# j is now the distance matrix for the graph.
for k in r:
    if j[k][-1]>n:
        # This means that vertex k is not connected to the "drain" vertex, and is therefore flooded.
        i[k/w][k-w*(k/w)]='*'
for r in e:print(''.join(i[r]))

The idea here is that I construct a directed graph where each grid cell has its own vertex (plus one additional "drain" vertex). There is an edge in the graph from each higher valued cell to its neighboring lower valued cells, plus there is an edge from all exterior cells to the "drain" vertex. I then use Floyd-Warshall to calculate which vertices are connected to the "drain" vertex; any vertices that are not connected will be flooded and are drawn with an asterisk.

I don't have much experience with condensing Python code, so there's probably a more succinct way I could have implemented this method.

Haskell, 258 characters

a§b|a<b='*'|1<3=a
z=[-1..1]
l m=zipWith(§)m$(iterate(b.q)$b(\_ _->'9'))!!(w*h)where
 w=length m`div`h
 h=length$lines m
 q d i=max$minimum[d!!(i+x+w*y)|x<-z,y<-z]
 b f=zipWith(\n->n`divMod`w¶f n)[0..]m
 (j,i)¶g|0<i&&i<w-2&&0<j&&j<h-1=g|1<3=id
main=interact l

Example run:

$> runhaskell 2638-Lakes2D.hs <<TEST
> 8888888888
> 5664303498
> 6485322898
> 5675373666
> 7875555787
> TEST
8888888888
566*****98
6*85***898
5675*7*666
7875555787

Passes all unit tests. No arbitrary limits on size.

• Edit (281 → 258): don't test for stability, just iterate to the upper bound; stop passing constant argument m

Python, 483 491 characters

a=dict()
def s(n,x,y,R):
 R.add((x,y))
 r=range(-1,2)
 m=set([(x+i,y+j)for i in r for j in r if(i,j)!=(0,0)and(x+i,y+j)not in R])
 z=m-set(a.keys())
 if len(z)>0:return 1
 else:return sum(s(n,k[0],k[1],R)for k in[d for d in m-z if a[(d[0],d[1])]<=n])
i=[list(x)for x in input().strip().split('\n')]
h=len(i)
w=len(i[0])
e=range(0,w)
j=range(0,h)
for c in[(x,y)for x in e for y in j]:a[c]=int(i[c[1]][c[0]])
for y in j:print(''.join([('*',str(a[(x,y)]))[s(a[(x,y)],x,y,set())>0] for x in e]))

I'm pretty sure there's a better (and shorter) way of doing this

Common Lisp, 833

(defun drains (terr dm a b)
  (cond
    ((= (aref dm a b) 1) t)
    ((= (aref dm a b) -1) nil)
    ((or (= a 0) (= b 0)
     (= a (1- (array-dimension terr 0)))
     (= b (1- (array-dimension terr 1)))) t)
    (t (loop for x from -1 to 1
       do (loop for y from 0 to 1
           do (if (and (or (> x 0) (> y 0))
                   (drains terr dm (+ a x) (+ b y))
                   (<= (aref terr (+ a x) (+ b y))
                   (aref terr a b)))
              (progn
                (setf (aref dm a b) 1)
                (return-from drains t)))))
    (setf (aref dm a b) -1)
    nil)))

(defun doit (terr)
  (let ((dm (make-array (array-dimensions terr))))
    (loop for x from 0 to (- (array-dimension terr 0) 1)
       do (loop for y from 0 to (- (array-dimension terr 1) 1)
         do (format t "~a"
            (if (drains terr dm x y)
                (aref terr x y)
                "*"))
         finally (format t "~%")))))

No attempt has been made to golf this, I just found the problem interesting. Input is the 2D array of the map. The solution checks each square to see if it "drains" -- a square drains if it is on the outer edge or if it is adjacent to an equal or lower height square that drains. To keep from recursing endlessly, the code keeps a "drain map" (dm) where it stores the drainage status of squares that have already been determined.