In geophysics determining Earth's thermal structure is key information allowing the development of evolutionary models and understanding the relationship between shallow and inner processes. The nature of the problem restricts obtaining samples and performing direct observations or tests to characterize the material. Therefore, this information is obtained by solving an inverse problem. Up to the depth of interest, the domain is composed of two layers: the \textit{Lithosphere} above that is rigid and colder, and the \textit{Astenosphere} below that is hot and flowing viscously. These layers are separated by the so-known Lithosphere-Asthenosphere Boundary (LAB), which is commonly assumed to be an isotherm. The immediate consequence of this definition is the presence of an essential condition immersed in the domain. Considering the existing boundary conditions, the problem results over-constrained. While this is well-known in the community, and several authors proposed different approaches to circumvent it, the strategies usually involve non-physical procedures, for example combining two independent temperature fields made compatible by smearing out the differences. Consequently, the solution in this portion does not comply with the governing equation. Considering the importance of inverse problems in geophysics, this work aims to develop a suitable tool for use in this context. With this aim, the essential condition introduced by the LAB isotherm is imposed weakly using Nitsche's method, ensuring it for both conforming and non-conforming meshes. The idea is to solve a pure diffusive heat flux problem in both sub-domains in successive stages. The strategy exploits the fact that boundary conditions are known with different degrees of confidence. Therefore, the idea is to relax the least precisely known boundary condition, such that the obtained solution, both in terms of temperature and flux, is as continuous as possible. The upper sub-domain boundary conditions are precisely known, and it is solved. Then, the difference between heat flux from the upper and lower sub-domain is minimized. Operationally requires expressing the lower sub-domain flux in terms of the unknown boundary condition to solve a minimization problem. As a result, stable and oscillation-free temperature fields are obtained. Moreover, the fluxes at the bottom of the domain are reasonable and compatible with the expected values.