By Jaan Kalda (TalTech).
For convenience, we assume T_0=0 as we can always use a temperature scale where T_0=0.
Since heat transport along the pipe is fast, and diffusion in the soil becomes increasingly slower over time, the radial temperature gradient dominates over the longitudinal one. Therefore, the heat flux density can be written as J=\kappa \frac{\partial T}{\partial s}, where s denotes the distance from the axis of the pipe. This means that the heat accumulated in an interval from s to s + \mathrm{d}s during \mathrm{d}t per unit length along the axis of the pipe equals \rho c (2\pi s \mathrm{d}s) \mathrm{d}T on one hand, and
\mathrm d(2\pi sJ)\mathrm dt=\kappa \frac{\partial}{\partial s}2\pi s\frac{\partial T}{\partial s}\mathrm dt\mathrm ds,
on the other hand. Once we equate these two expressions, we obtain
\frac{\partial T}{\partial t}=\frac{\kappa}{\rho cs} \frac{\partial }{\partial s}s\frac{\partial T}{\partial s}.\;\;\;\;\;\;\;\; (1)
This is a diffusion equation in cylindrically symmetric geometry, with the ratio D = \frac{\kappa}{\rho_w c_w} referred to as the diffusivity. As discussed by the hints, the characteristic time for heat to propagate over a distance \delta is estimated as \delta^2/D (this estimate can also be obtained from a dimensional analysis). Hence, we can expect that if the system evolution time is larger than \delta^2/D, a quasi-equilibrium temperature distribution is reached over regions of size \delta; equivalently, if the waiting time is t, then a quasi-equilibrium is reached in regions with size smaller than \sqrt{4Dt} (the factor 4 is here due to the fact that the root-mean-square displacement of Brownian particles in 2D geometry is \sqrt{4Dt}). This means that the left-hand side of the diffusion equation is approximately zero: \frac{\kappa}{\rho c s} \frac{\partial}{\partial s} \left( s \frac{\partial T}{\partial s} \right) = 0. Hence, the heat flux per unit length of the pipe P = -2\pi\kappa s \frac{\partial T}{\partial s} is constant, and we obtain \frac{\partial T}{\partial s} = -\frac{P}{2\pi\kappa s}, hence
T=-\frac P{2\pi\kappa}\ln \frac sr+T_\star,
assuming that s\ll\sqrt{4Dt}. Using the condition
T(s=r)=T_1
we obtain the value of the integration constant T_\star=T_1.
The constant P can be found from the condition that the heat flux P supplies heat necessary to expand the heated region around the pipe. When time increases from t to t + \mathrm{d}t, the increase in the cross-sectional area of the heated region is on the order of 4\pi D \mathrm{d}t, and the peripheral region is heated up to the temperature on the order of T(s=\sqrt{4Dt}). Therefore, the heat balance can be written as:
P\mathrm dt\approx 4c\rho\pi D\mathrm dt(T_1-\frac P{2\pi\kappa}\ln \frac {\sqrt{4Dt}}r)=4\pi\kappa \mathrm dt T_1-P\mathrm dt\ln \frac {4Dt}{r^2}.
From here, we can conclude that
P\approx\frac{4\pi\kappa T_1}{\ln(\mathrm 4\mathrm eDt/r^2)}.\;\;\;\;\;\;\;\; (2)
Note that since we are only estimating the order of magnitude, minor differences in estimates (e.g. dropping the factor 4, or not having the factor e) are OK.
Now, with this result at hand, we can easily finalize the solution by equating the heat loss to the soil over the entire length of the pipe, PLt, to the heat supplied by hot water:
tL\frac{4\pi\kappa T_1}{\ln(4\mathrm e\kappa t/\rho cr^2)}\approx \pi r^2uc_w\rho_wT_1t,.
Hence,
t\approx \frac{r^2c\rho}{4\mathrm e \kappa}\exp\left(\frac{4L\kappa }{r^2uc_w\rho_w}\right).\;\;\;\;\;\;\;\; (3)
Pay attention to the fact that at the last step, we assumed that the evolution time t in Eq. (2) is almost constant over the entire length of the pipe, so that the total loss into the soil can be obtained by multiplying P with the length L. Strictly speaking, the time t during which the heat has propagated radially away from the pipe decreases with the increasing distance from the hot water inlet. However, we know now that the distance of the temperature front from the inlet grows logarithmically in time. In other words, the front speed becomes slower and slower, and most of the time is spent waiting for the front to pass the very last segment of the total length L. This means that the evolution time t in Eq. (2) is, indeed, almost constant over the entire length of the pipe.
More precise calculations
The accuracy of the estimate (2) can be improved by obtaining a self-similar solution of the diffusion equation (1) which satisfies approximately the boundary conditions. Finding self-similar solutions is a standard technique when dealing with equations in partial derivatives (PDEs). A self-similar solution is a solution that can be rescaled at any moment in time to its original shape. In our case, we look for the solution of Eq. (1) in the form
T(s,t)=s^\alpha T(s/t^\beta).\;\;\;\;\;\;\;\; (4)
The value of the exponent \alpha is determined by the initial and boundary conditions. Our boundary condition T(s=r, t)=T_1 cannot be handled by such a self-similar solution, but let us keep in mind that r is small. If we were able to set r=0, the obtained boundary condition T(s=0,t)=T_1 would correspond to \alpha=0. Strictly speaking, we cannot set r=0 because, as it appears, the solution will diverge logarithmically at s=0. However, logarithmic divergence is slow, and for the time being, let us accept \alpha=0 without worrying about the divergence at s=0.
The value of the exponent \beta is determined by the structure of the PDE: it needs to be such that if we substitute (4) into (1), the resulting ordinary differential equation will be such that t and s enter only in combinations substitutable with \xi=s/t^\beta, i.e., the obtained ordinary differential equation is an equation for T(\xi). We can easily see that we need to take \beta=\frac{1}{2}. If we take a dimensionless variable
\xi=s/\sqrt{2Dt},
we will be able to get rid of all the factors in our equation for T(s,t)\equiv T(\xi):
(\xi T')'+\xi^2T'=0.
If we introduce a new function y(\xi)=\xi T', our equation simplifies to
\mathrm dy+y\xi\mathrm d\xi=0
which can be easily solved: y = a \exp(-\xi^2/2). This means that T = \int y \mathrm{d}\xi/\xi = a \int \exp(-\xi^2/x) \frac{\mathrm{d}\xi}{\xi} = \frac{a}{2} \mathrm{Ei}(-\xi^2/2) + b, where \mathrm{Ei}(x) is the exponential integral, and b is an integration constant. Since \mathrm{Ei}(-x) tends to zero as x \to \infty and the same needs to hold for T(x,t), we can conclude that b = 0. Therefore, T = \frac{a}{2} \mathrm{Ei}(-r^2/4Dt), where the factor a can be determined from the condition T(s=r) = T_1. To sum up,
T=T_1\frac {\mathrm{Ei}(-r^2/4Dt)}{\mathrm{Ei}(-s^2/4Dt)}..\;\;\;\;\;\;\;\; (5)
There is one problem with this solution: it is not precise, because the factor a needs to be constant, but in our case, it is not. To make it precise, we can freeze the prefactor:
T=T_1\frac {\mathrm{Ei}(-r^2/4Dt)}{\mathrm{Ei}(-s^2/4Dt_0)}.\;\;\;\;\;\;\;\; (6)
where t = t_0 is the current moment of time, but in that case, our solution satisfies a slightly different boundary condition
T(r=s)=T_1\frac {\mathrm{Ei}(-s^2/4Dt)}{\mathrm{Ei}(-s^2/4Dt_0)}.
Here the argument is very small, s^2/4Dt \ll 1, so we can use a series expansion for the exponential integral,
\mathrm{Ei}(-r^2/4Dt)\approx\gamma-\ln(4Dt/r^2),
where \gamma \approx 0.57721 is the Euler–Mascheroni constant. Thus, the boundary condition can be rewritten as
T(r=s)=T_1\frac {\ln(4Dt/s^2)-\gamma}{\ln(4Dt_0/s^2)-\gamma}.
Since the asymptotic growth rate of the logarithm is very slow, the right-hand side of this expression becomes almost constant for t, t_0 \gg s^2/4D, i.e., near t = t_0, Eq. (6) almost satisfies the boundary condition, and we can expect that for t = t_0 \to \infty, the true solution of Eq. (1) is asymptotically approaching Eq. (5), which can be approximated for r \ll 2 \sqrt{Dt} as
T=T_1\frac {\ln(4Dt/r^2)-\gamma}{\mathrm{Ei}(-s^2/4Dt)},
and for r,s\ll 2\sqrt{Dt} as
T=T_1\frac {\ln(4Dt/r^2)-\gamma}{\ln(4Dt/s^2)-\gamma}.
Now we are ready to work out an improved version of expression (2) for the heat loss per unit length of the pipe:
P=2\pi\kappa r\left.\frac{\mathrm dT}{\mathrm ds}\right|_{s=r}=\frac {4\pi\kappa T_1}{\ln(4Dt/r^2)-\gamma}.
Next, we write down the heat balance equation along the pipe,
\mathrm d(\pi r^2c_w\rho_wvT)=-P\mathrm dz=-\frac {4\pi\kappa T}{\ln(4Dt/r^2)-\gamma}\mathrm dz,
to obtain the water temperature as a function of the coordinate z along the pipe:
T=T_1\mathrm e^{-z/\Lambda},\;\; \Lambda=[\ln(4Dt/r^2)-\gamma]\frac{ r^2c_w\rho_wv}{4\kappa}.
We need to have T(z=L) = T_1/2, hence L = \Lambda \ln 2 which yields
\ln(4Dt/r^2)=\gamma+\frac{4\kappa L}{ r^2c_w\rho_wv\ln 2},
and finally
t=\frac{r^2c\rho}{4\kappa}\exp\left(\gamma+\frac{4\kappa L}{ r^2c_w\rho_wv\ln 2}\right).