What is the physical interpretation of the Wick rotation? What is the physical interpretation of the Wick rotation?
How is it that we can just propose there's a new time coordinate tau? Are physicists saying time is modeled by an imaginary number? Isn't that at odds with the clocks we use everyday?
 A: There's a couple of different points of view on what exactly it means.
Wick Rotation and diffusion
Even in the Schrodinger equation, Wick rotation has an interesting viewpoint. Namely, that Wick rotating the Schrodinger equation gives us the diffusion equation! More explicitly, the Schrodinger equation (say, with no potential) is given as
$$i\partial_t \psi + \nabla^2\psi = 0$$
Wick rotation of $t = i\tau$ gives 
$$\partial_t \psi - \nabla^2\psi = 0$$
which is precisely the diffusion equation. So quantum propagation in real time is, roughly-speaking, diffusion in imaginary time! 
Boltzmann statistics and the Hamiltonian
Remember from statistical mechanics that the partition function of a system at inverse temperature $\beta = 1/{k_B T}$ is given by
$$Z(\beta) = \Sigma_i e^{-\beta E_i}$$.
It turns out that this thermal quantity can be defined in terms of the Hamiltonian, namely that 
$$Z(\beta) = Tr(e^{-\beta H})$$
where $\rho(\beta) = e^{\beta H}$ is the matrix known as the density matrix of the Hamiltonian. This density matrix is actually all we need to provide a complete description of the thermal properties of a quantum system. 
Also note that $\rho(\beta)$ is closely tied together the time-evolution operator $U(t) = e^{iHt}$ by setting $\beta = -it$. This is closely tied together with what we said about diffusion earlier. Remember that a thermal system is an equilibrium notion about letting a system relax given some energy constraints. But the abstract mathematical standpoint is telling us that letting a system thermalize to an inverse temperature $\beta$ is the same thing as letting it evolve in imaginary time for a time $\beta$. But from before, we can see that evolving in imaginary time is the same as letting the system diffuse for a time $\beta$.
In particular, we can also say something about the ground state of the quantum system from this viewpoint. Mathematically speaking, if there is a unique ground state $|\psi>_{ground}$ at energy zero, then it can be obtained (in an unnormalized form) from any initial state $|\psi>_{0}$  with nonzero overlap from the ground state:
$$|\psi>_{ground} = \lim_{\beta \rightarrow \infty} e^{-\beta H} |\psi>_{0} $$
(which is a simple exercise that one should do.) So, from this viewpoint, the ground state of a system is obtained by letting the system evolve for an infinite amount of imaginary time from any starting point, or by letting it "diffuse" for an infinite amount of time. 
Note that I've necessarily been a bit hand-wavy here. No one truly knows what the precise connection between diffusion and real-time evolution is apart from the mathematics, but presumably, it will be more clear once we understand what quantum mechanics truly is. 
From a computational standpoint
Wick rotation just transforms a path integral into something that actually converges. To understand this point, it suffices to just look at a simple one dimensional integral
$$I(\alpha) := \int_{-\infty}^\infty e^{i \alpha x^2} dx$$
For $\alpha$ being a purely a real number, this integral does not converge, since the integrand has norm 1 at all $x$. But, if we shift $\alpha \rightarrow \alpha + i \epsilon$ for some small number $\epsilon$, then we the imaginary part forces the integrand to go to zero exponentially as $x\rightarrow \infty$, and will converge to $I(\alpha) = \sqrt{\frac{\pi}{-i \alpha}}$. From here, we see that we can evaluate the non-convergent integral $I(\alpha)$ by simply making $\alpha$ into a slighly imaginary number, which allows us to define the integral for $\alpha$ purely real.
The reason that this is important for quantum field theory is that the path integral 
$$\int D\phi e^{iS[\phi]}$$
is an infinite dimensional integral. Moreover, a prototypical example of such an action looks like
$$S[\phi] = \int d^4x \phi(x)(-\partial_t^2 + \partial_x^2 + m^2)\phi(x)$$
The comparison point between this infinite dimensional integral and the finite dimensional one is that they are both quadratic matrices with imaginary eigenvalues all with real part bigger than zero. (The spectrum of the operator $(-\partial_t^2 + \partial_x^2 + m^2)$ turns out to satisfy this property.) So, evaluating the path integral in this case will give us an infinite number of nonconvergent integrals! Luckily, making the transformation $t = -i\tau$ will make the integral in terms of $\tau$ a convergent integral for real $\tau$ the lesson learned from the 1D case will allow us to evaluate it for a small imaginary.
The reason this silly-looking substitution works is because QFT's are holomorphic. This is typically an axiom of QFT's and means that any observable $<{O(x_1)...O(x_n)}>$ is a holomorphic function of the $x_i$. So, making time imaginary and continuing back to real time allows us to actually define finite quantities in the first place. An important caveat of this is that there are generically poles in these correlation functions, so we actually have to be careful in how we analytically continue.
A: Besides the nice answer by Joe, I would like to mention a noncommutative-geometry-inspired approach -- which is more 'geometric' in nature than 'physical' and attempts to arrive at an interpretation of the Wick rotation. Disclaimer: Now, one could, in principle, argue that 'geometry is physics' as in the sense of Einstein: 'geometry is gravity'. So take how 'physical' it is, at all, with a pinch of salt, as per your taste.
Noncommutative geometry à la Connes is defined by a spectral triple, which could be twisted by algebra automorphisms [1], [2]. Although the motivation for the twist was mathematical: to construct spectral triples for type III algebras [3], it quickly found applications in quantum groups, C*-dynamical systems, Standard Model and its extensions, cf. [4] and references therein.
The twisted noncommutative Standard model was intended for electroweak vacuum stability and fixing Higgs mass. Later twisting was found to naturally yield a Krein space (of Lorentzian spinors) and hinted at a tight link with Wick rotation [5]. In a recent paper, twisting Euclidean spectral triples and evaluating fermionic action induced on them by the twist has led to the derivation of Weyl and Dirac equations in Lorentzian signature [6].
In this setting, Wick rotation is something that naturally appears, at least at the action level, as a consequence of twisting Euclidean spectral triples. One might even go ahead and say that the transition to a Lorentzian geometry from a Euclidean one, in this sense, is a result of a (spontaneous?) symmetry-breaking: from a non-twisted geometry to a twisted one.
