Imagine digging a cylinder-shaped (vertical) bore-well of depth $L$ and diameter $r$ ($L\gg r$). The (infinitely thin) cylinder-wall is made watertight and we split the well in half using a kind of thin lamella through the central axis of the cylinder. On one side of this lamella we pump water down, allowing it to come back up on the other side. The goal of this operation is to get some useful geothermal heat out of the crust into this water. Imagine for simplicity that our goal today is to maximize the amount of heat we want to extract in that operation (as opposed to e.g. getting a high temperature in the returning water) and therefore we will push the water pumping rate to the limit where the water doesn't get time to heat up significantly during such a cycle and so we assume the borewell is at a constant temperature $T=0$ (the surface temperature). Assume that the crust conducts heat according to the stationary heat equation (Laplace eq.) with constant heat conductivity (energy current $J=-\nabla T$ which means the temperature profile in the crust, before we have dug our well, will be linear along the vertical direction in this model: $T_0(x,y,z)=z$ understanding that $z\uparrow$ means going deeper and consistent with the temperature gauge which puts the temperature at the surface at zero) and assume that after digging the borewell and pumping the water the surface temperature will still remain $0$. What is the (stationary) rate at which we could extract heat from the crust into this bore-well? (The question is meant to gather more elegant, efficient answers than what my own thinking has produced so far)
I went about estimating this quantity by superposing onto the original temperature profile $T_0$ the Coulomb potential* associated to the line charge $\rho(x,y,z)=-\alpha\delta(x)\delta(y)\chi_{[-L+r,L-r]}(z)z$ where $\alpha=[2\sinh^{-1}(L/R-1)-2]^{-1}$ (supposing that the borewell's central axis is the $z$-axis) and this led to the following lower (but probably relatively sharp) bound for the absorbed heat: $$\frac{\pi (L-r)^2}{\sinh^{-1}(L/r-1)-1}(=-4\pi\int_{\mathbb{R}^2}dx\,dy\int_{\mathbb{R}^+}dz\,\rho(x,y,z))$$ The reason why it's a lower bound is that the final temperature profile resulting from this 2-step construction has a boundary temperature greater than $0$ on the well while the exact temperature profile would have boundary value $0$ there. By subtracting the former from the latter temperature profile, observing that this difference solves the Laplace equation it is possible to infer that the former T-profile carries less heat into the bore-well than the latter.
*With "Coulomb potential $\Delta T$ associated to $\rho$" I mean $$\Delta T(x,y,z)=\int_{\mathbb{R}^3}dx'dy'dz'\frac{\rho(x',y',z')}{((x-x')^2+(y-y')^2+(z-z')^2)^{1/2}}$$
