Non-convergance when calculating temperature/heat flows through a section of rock
I am attempting to calculate temperature of section of rock in the earth as a function of vertical position in the rock and time. Along with it I am calculating the heat flow through the rock as a function of vertical position in the rock and time. Note that this is a one-dimensional examination of the problem- imagining that the rock is uniform in every direction but up-down.
Please be patient with me: I currently do not have the vocabulary nor the background to talk at the level of familiarity with FEM and differential equations I see in other questions on this site. I am a software developer with some physics background, which has allowed me to get this far.
Description of methodology: I am currently splitting up the section of rock into horizontal slices (depth steps), and then calculating the temperature and heat flow of each depth step at specific time-steps. I am assuming that heat only flows upward from below the section of rock.
As I understand it, the temperature difference between the top and bottom of a depth step would be expressed as: $$ \Delta T = \frac{Q\Delta Z}{\lambda} $$ Where $Q$ is the heat flow at the middle of the depth step in $mW/m^2$, $\Delta Z$ is the thickness of a depth step at this time, and $\lambda$ is the thermal conductivity of the rock (constant).
Additionally, the change in the heat flow between the bottom and top of a depth step could be expressed as: $$ \Delta Q = - C_p \rho \Delta Z \frac{\Delta T}{\Delta t} $$
where $C_p$ is the heat capacity of the rock (constant), $\rho$ is the density (constant), and $\frac{\Delta T}{\Delta t}$ would be the change in temperature of the depth step over the course of this time-step.
As a starting point for the calculations, I know the temperature at the top of the rock layer as a whole (but not the temperature at the bottom of the rock layer). I also know the heat flow value just below the rock layer (but not the heat flow value at the top).
The key piece of complexity to this is that the rock layer is increasing in thickness through time (aka deposition), so that I am adding to the $\Delta Z$ of the top-most depth step and then adding more depth steps to the top of the section of rock as time progresses. As a consequence, individual depth steps get more and more buried, and thus hotter through time, which thus causes the heat flow to decrease as you move upwards in the rock layer (as per equation 2 above).
My method of calculating the temperatures and heat flows for each time-step is as follows:
1. Increase the thickness of the rock layer, adding depth steps to the top as necessary.
2. 749
Calculate the temperatures at the top and bottom of each depth step going down the layer (since I know the surface temperature, and thus the starting point), using the heat flow values from the previous time-step.
3. Calculate the heat flow values at the bottom and top of each depth step going up the layer (since I only know the heat flow at the very bottom of the rock layer), using the temperatures of each depth step calculated in step #2, along with the temperature in the previous time-step.
4. Repeat steps 2-3 until the solution converges. Note that in step #2 I can now use the heat flows calculated in step #3.
The Problem: This solution does not always converge. As the rock layer gets thicker and thicker through time, each time-step takes longer and longer to converge, until it reaches a tipping point and starts to diverge. What can I do to force this to converge?
I am suspicious that the reason for divergence is because I'm forced to calculate temperature down the depth steps and heat flows up the depth steps. This means that the bottom-most depth step's temperature is sensitive to ALL the changes in temperature above it as per the first equation above. But this then causes the heat flow at the bottom to change from the previous iteration by a fair amount, which then propagates all the way to the top of the rock layer via the second equation above, which then affects the temperature at the top strongly, etc.
Warning: The above problem is actually a major simplification of a larger problem I'm trying to solve- the constants in the equation above are actually temperature dependent and depth dependent in a complicated, non-continuous way. However, the problem I am facing still exists in the simplification as described above. 108