It sounds simple but requires at least some basic knowledge of partial differential equations. You also need to know about discretization so that you can solve the equation
\[ cp \rhow\frac{\partial T}{\partial t} = K \frac{\partial^2 T}{\partial x^2} \]
on your computer. Most commonly used are finite difference schemes, for example, explicit (forward time stepping), implicit (backward time stepping), or the Crank-Nicolson (both, forward and backward time stepping) schemes.
Update:
Here you see an example of a vertical temperature profile of a snowpack of 5m thickness. On top (left, 0m) the model is forced with 10 years of surface temperature at location of the JAR automatic weather station as modelled by MAR. You can see the alternating summer and winter heat and cold waves entering the snowpack. The annual cycle is the strongest because a slower surface signal (i.e., longer time scales) lead to a larger penetration depth \(l\) \[ l = \sqrt{2\,\frac{K}{cp \rhow} \Delta t} \]
Because of the (almost) periodic forcing the signal at depth is phase shifted and decays exponentially.
And I like coding. I've been working with complex computer models ever since my undergrad and I enjoy data exploration and data analysis to gain insights into the underlying principles.
Feel free to contact me.