# How to implement boundary conditions in heat equation with no flux and fixed value at the same time? Is it Robyn BC?

I am modeling the temperature of the groundwater using heat equation. I have Dirichlet BC at the top but at the bottom I have constant temperature equal 12 degrees C (see attached pic). It is look like at the bottom we have zero flux. When I specify the Dirichlet BC I have a "knee" right before "neutral zone" which is not the case. So far I have 2 ideas how to implement it:

1. I was thinking maybe the Robyn BC will help here?
2. Maybe I need to implement somehow the [-hh * (T - T_inf)/k] as gradient there. When we have not 12C we will have the flux and when we have 12 then we will have 0 gradient.

what do you think? is it correct way of tackling this problem?

• This stackexchange isn't really the right place to ask about the correct mathematical statement of the problem. You need to know that before you try computing anything. – David Ketcheson Oct 7 '15 at 7:36
• I'm voting to close this question as off-topic because it is about the physics or engineering of the problem, not about its computational solution. – David Ketcheson Oct 7 '15 at 7:37

You cannot (in general) both enforce zero flux and constant value at the same boundary.

Additionally, you don't need a Robin boundary condition. The boundary condition depends on what you wish to enforce. At the top, a Dirichlet condition seems natural. You have already suggested two options at the bottom: Dirichlet (constant, 12 degrees as you said) or Neumann (zero flux, as you said).

In this problem, if the bottom boundary condition is a great enough distance from the top there will be only a very small difference between a zero flux Neumann condition and a Dirichlet condition at the mean value of the annual surface temperature. How far is "a great enough distance" depends on the skin depth as mentioned by @GeoMatt22. If you are seeing a bend just above the "neutral zone" with the Dirichlet condition it is probably best to switch to Neumann.

Do you have variable thermal diffusivity with depth?

If not, then if I understand you problem correctly, there is an analytical solution (see here).

In either case, assuming I understand your problem correctly, then typically the lower boundary condition is zero heat-flux "at infinity" (i.e. Neumann). Here "infinity" means at $\frac{z_0-z}{\lambda}\gg1$, where $\lambda$ is the "skin depth" (see the above link). With this condition, the temperature will end up as the time-average of your surface-forcing.

Edit: Your comment suggests you have a piecewise-constant diffusivity, dry soil vs. saturated soil (ignoring the vadose zone). Across a step change in diffusivity, the temperature $T$ and heat flux $q=-KT_z$ will be continuous, but the temperature gradient $Tz$ will be discontinuous. Is the kink in your results at the interface (i.e. water table)?

Numerically, piecewise constant conductivities would be handled in a finite volume method by assigning the interface-conductivity as the harmonic average of the two adjacent cell-conductivities when computing the flux. (These notes give an overview.)

• yes, I have. I have water at the bottom and air at the top of the soil. – Igor Markelov Oct 6 '15 at 14:00