Does Riemann map depend continuously on the domain?

Here's a conceptual proof why this is true, up to things which are intuitively obvious and not hard to prove:

In the unit disk, almost every Brownian path hits the boundary. The hitting measure equals proportional to arc length. In two dimensions, a conformal map takes trajectories of Brownian paths to trajectories of Brownian paths: just the time parametrization changes. (This is a consequence of the fact that conformal maps take harmonic functions to harmonic functions; harmonic functions are the functions whose expectation is invariant under Brownian motion.)

It follows that the pushforward of arc length of the unit disk via the Riemann mapping is the hitting probability for Brownian paths starting at the image of the origin. Your question is equivalent to asking whether the measure of intervals in your parametrized Jordan curves is uniformaly continuous with respect to the uniform topology on parametrized Jordan curves.

It's intuitively obvious, as well as true and not hard to prove (further explanation below), that a Brownian path that starting at a point $z$ inside a Jordan domain near the boundary is likely to hit the boundary nearby. This fact quickly implies the continuity that you need: follow Brownian motion until it gets within $2 \epsilon$ of the initial boundary curve, so it is between $\epsilon$ and $3 \epsilon$ of the perturbed curve. When Brownian motion is continued, most of it can't shift very far before hitting.

(Note: given a Jordan curve, you must take $\epsilon$ small enough that short intervals as measured by hitting measure are also short on the curve, to be able to conclude that the Riemann mapping does not move very far when you perturb the curve.)

There are a variety of ways to prove that a ; one way is to lift continuously to a branch of the map $\log(z-z_0)$, where $z_0$ is the closest point to the boundary. Now the random walk takes place in an arbitrarily long strip of width no more than $2 \pi$; it has little chance of remaining in the strip long enough to move far along its length. This follows from the fact that a Brownian path in 1 dimension has a large probability of going outside an interval of length $2 \pi$ after a certain length of time.

Another way to prove that Brownian paths are likely to hit nearby on the curve is to make use of the estimate for the Poincaré metric inside a domain: it varies by no more than a factor of 2 from 1/(minimum distance to the boundary). With this estimate, you can show that for a large Poincaré disk centered about $z$ near a boundary point $z_0$, most of its arc length gets squeezed near to $z_0$.

Side note: Brouwer proved (in his intuitionistic framework) that every function that is everywhere defined is continuous, so from this point of view Caratheodory's theorem about continuity at the boundary implies continuity. However, one needs to check that Caratheodory's theorem is true intuitionistically; Brouwer later rejected his famous fixed pointed theorem on these grounds.

Yes. This is a classical result of the geometric function theory due to Carathéodory.

The theorem and a fairly straightforward proof can be found, for instance, in the Hurwitz-Courant Funktionentheorie (in the part written by Courant).

Edit. Theory of Functions of a Complex Variable by Markushevich might be a more accessible reference (see Volume 3, Theorem 2.1).

Edit 2. By the way there exist sequences of domains $\{G_n\subset\mathbb C:\ n\in\mathbb N\}$ such that $$\limsup_{n\to\infty}\ \mbox{dist}(\partial G_n, \partial G)>0$$ but the corresponding Riemann maps $\phi_n:\mathbb D\to G_n$ converge uniformly in $\mathbb D$ to the Riemann map $\phi:\mathbb D\to G$.

For example, let $G_n$ be a union of two disjoint rectangles $G'$ and $G''$ connected with a 'thin' rectangle of the fixed length $l$ and width $h_n=1/n$ (see the picture below).

Let $z_0\in G'$. Then the sequence of the conformal maps $f_n:G_n\to\mathbb D$ which satisfy the conditions $$f_n(z_0)=0,\qquad f'_n(z_0)>0,$$ converges uniformly in $G'$ to the conformal map $f:G'\to\mathbb D$ satisfying the same condition. Moreover, the sequence of the inverse maps $\phi_n:\mathbb D\to G_n$ converges uniformly in $\mathbb D$ to $\phi=f^{-1}$. The general Carathéodory theorem gives a criterion of the convergence of the Riemann maps in terms of the corresponding domains.

This fact is precisely Theorem 2.11 in the book Boundary Behaviour of Conformal Maps by Pommerenke. It follows from Caratheodory's Kernel Theorem.