Closed plovely closed 5 years ago
The model has no radiogenic heat production, so in steady-state the heat flux should be constant everywhere below the land-air interface.
The model has no radiogenic heat production and the diffusivity is constant, at steady state, the gradient of the geotherm should be constant...
In this model, diffusivity is homogeneous, as is heat capacity, so thermal conductivity should be proportional to density (k=alpha rho Cp).
Yes I agree with that...
Because density of the mantle is greater than density of the crust, conductivity of the mantle should be greater and, given q=-k * grad(T), thermal gradient should be less in the mantle than in the crust. Instead, what I see (figure below) is a constant gradient in the temperature field and relatively greater heat flux in the mantle than in the crust.
Because density is higher, conductivity of the mantle should be higher (assuming constant diffusivity) and as the gradient is constant, the heat flow should be higher in the mantle than in the crust... I think what you see is correct.
OK, maybe I'm missing something obvious here with regard to thermal conductivity vs. diffusivity.
I'm looking at equation 4-45 in Turcotte & Schubert (2nd edition, pg 146) for steady state heat conduction in 2D: dq_x/dx + dq_y/dy=rho*H In this case, H=0 (no radiogenic heat production), and it's effectively a 1-D problem so dq_x/dx=0. However, the discontinuity in heat flux from mantle to crust means that dq_y/dy is not equal to zero everywhere. In other words, thermal energy is accumulating at the discontinuity in q and temperature isn't in equilibrium
Romain, Julian,
I think that this issue needs to be reopened. Upon comparison of UW models with other thermal models, I'm observing non-trivial discrepancies in pre-rift (steady-state) heat flow.
In response to Romain's comment from 24 April, my issue is that heat flow in (at the base) is not equal to heat flow out (at the earth surface). Therefore, in the absence of an internal heat source, the model is not in steady state.
To conserve energy in steady state, the geothermal gradient should be constant across layers with equivalent conductivity, not necessarily across layers with equivalent diffusivity. I have run another pseudo 1-D model in which heat capacity is constant, and the product of diffusivity and density are constant. Thus, conductivity is homogeneous. However, thermal gradient and heat flow are discontinuous, and therefore I don't think the model is in steady state at the start:
I will share the notebook in a separate email.
Thanks, Peter
Hi Peter
Since the product diffusivity x density is constant in your model, all other parameters being constant between crust and mantle, then the heat flow discontinuity you observe doesn’t make sense.
Perhaps a default value is overwritting the inputed value?
Cheers
Patrice
On 10 Jan 2020, at 6:19 am, plovely notifications@github.com<mailto:notifications@github.com> wrote:
Romain, Julian,
I think that this issue needs to be reopened. Upon comparison of UW models with other thermal models, I'm observing non-trivial discrepancies in pre-rift (steady-state) heat flow.
In response to Romain's comment from 24 April, my issue is that heat flow in (at the base) is not equal to heat flow out (at the earth surface). Therefore, in the absence of an internal heat source, the model is not in steady state.
To conserve energy in steady state, the geothermal gradient should be constant across layers with equivalent conductivity, not necessarily across layers with equivalent diffusivity. I have run another pseudo 1-D model in which heat capacity is constant, and the product of diffusivity and density are constant. Thus, conductivity is homogeneous. However, thermal gradient and heat flow are discontinuous, and therefore I don't think the model is in steady state at the start: [image]https://user-images.githubusercontent.com/30320218/72097695-91e62c00-32e2-11ea-9760-b42a152d7a41.png
I will share the notebook in a separate email.
Thanks, Peter
— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHubhttps://github.com/underworldcode/UWGeodynamics/issues/118?email_source=notifications&email_token=ADI3DZJKFNVSQROFUDKVOADQ452DPA5CNFSM4HG62SB2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEIRONYI#issuecomment-572712673, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ADI3DZMCBYVSUFWEBVOU5NLQ452DPANCNFSM4HG62SBQ.
Thanks, Patrice. I reran the model without a default (global) diffusivity value, and results are the same. If a default is overriding the prescribed parameters, it must be in the background.
I'll forward the notebook and ParaView state files that I sent to Romain and Julian, if you wish to take a look.
Cheers, Peter
Hi Peter,
How did you make that graph? Are you using the initial densities and the diffusivity? Or is it calculated using the result of the Steady State solve. I am a bit confused here.
Romain
On 1/10/20 10:51 AM, plovely wrote:
Thanks, Patrice. I reran the model without a default (global) diffusivity value, and results are the same. If a default is overriding the prescribed parameters, it must be in the background.
I'll forward the notebook and ParaView state files that I sent to Romain and Julian, if you wish to take a look.
Cheers, Peter
— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub https://github.com/underworldcode/UWGeodynamics/issues/118?email_source=notifications&email_token=AA4CVFTOR74BRZZYQ4FCTSDQ46Z63A5CNFSM4HG62SB2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEISGFHA#issuecomment-572809884, or unsubscribe https://github.com/notifications/unsubscribe-auth/AA4CVFUQDWR4MHSVN2TZIFDQ46Z63ANCNFSM4HG62SBQ.
Hi Romain,
Sorry, it's not straightforward. Here are the steps I followed in ParaView. I would send a screenshot or two, but am having trouble logging in from home.
1) load fields.xdmf
2) calculate gradient of temp field
3) run calculators on the gradients:
a) HF_crust = -gradient/1000 projDensityField alpha_crust cp
b) HF_mantle = -gradient/1000 projDensityField alpha_mantle cp
i divide by 1000 to convert the gradient from K/km to K/m
alpha is diffusivity. alpha_crust is 7.5e-7 and alpha_mantle is 6.31.
Density for the crust is 2778 and density for the mantle is 3300
cp is 1200
4) use the Plot Over Line filter to create a vertical profile through the center of the model. Display HF_crust and HF_mantle.
5) In the image I sent you I blocked out the depth range of each curve that is irrelevant (i.e. where HF_crust is calculated in the mantle and where HF_mantle is calculated in the crust).
I hope this helps. If you're able to open the ParaView state file with the output of the notebook I sent, it should help, as it contains steps 1-4. If you need further description I'll follow up in the morning when I'm back to the office.
Cheers, Pete
Also, the plot is shown for the steady state solution (t=0).
Cheers, Pete
On Thu, Jan 9, 2020, 6:34 PM Romain Beucher notifications@github.com wrote:
Hi Peter,
How did you make that graph? Are you using the initial densities and the diffusivity? Or is it calculated using the result of the Steady State solve. I am a bit confused here.
Romain
On 1/10/20 10:51 AM, plovely wrote:
Thanks, Patrice. I reran the model without a default (global) diffusivity value, and results are the same. If a default is overriding the prescribed parameters, it must be in the background.
I'll forward the notebook and ParaView state files that I sent to Romain and Julian, if you wish to take a look.
Cheers, Peter
— You are receiving this because you modified the open/close state. Reply to this email directly, view it on GitHub < https://github.com/underworldcode/UWGeodynamics/issues/118?email_source=notifications&email_token=AA4CVFTOR74BRZZYQ4FCTSDQ46Z63A5CNFSM4HG62SB2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEISGFHA#issuecomment-572809884>,
or unsubscribe < https://github.com/notifications/unsubscribe-auth/AA4CVFUQDWR4MHSVN2TZIFDQ46Z63ANCNFSM4HG62SBQ .
— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/underworldcode/UWGeodynamics/issues/118?email_source=notifications&email_token=AHHKMWQT66NSIL4GPCO3M5LQ467CDA5CNFSM4HG62SB2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEISIW6A#issuecomment-572820344, or unsubscribe https://github.com/notifications/unsubscribe-auth/AHHKMWTCHS6GCM2HZEHRJWLQ467CDANCNFSM4HG62SBQ .
Hi All,
I hope that the issue I'm seeing is becoming clear, and that the root cause is in how I built the model or how I'm analyzing the results, not the UW solution itself.
I realized there is an easier way to create effectively the same plot as I shared last week... Simply multiply the temperature gradient (of the steady state solution at t=0) by thermal conductivity, which is uniform (2.5 W/m/K) for crust and mantle. Then extract the gradient along a vertical line through the center of the model.
Here is a screenshot from ParaView. Note in the calculator that factors of 1000 to convert km to m (in the gradient) and to convert W to mW cancel.
Cheers, Peter
Hi Peter, Happy new 2020!
The steady state solve can be interpreted as zeroing the temperature field variation in time, ie, dT/dt = 0. The equation we use assumes velocity is 0 and is expressed with diffusivity, K, on the L.H.S, ie.
grad( Kgrad(T) ) = H/(rhocp)
So with H = 0 variations in rho don't affect the steady state equation because the R.H.S. = 0. Only variations of material diffusion will lead to a change in gradient of the geotherm.
In the model you sent via email the thermal gradients are discontinuous because the diffusivity values are also discontinuous, mantle.diffusivity != uppercrust.diffusivity.
Hope this helps.
Hi Peter,
I’d also suggest to make such plots without paraview. Paraview is great for quick visualisation and analysis but it could be a bit of a simplified plotting method if you’re interested in the true values of the function you want to plot. Matplolib or plotly could be reasonable alternatives you can use for in-line or post processing analysis.
Sent from my iPhone
On 15 Jan 2020, at 9:20 am, Julian Giordani notifications@github.com wrote:
Hi Peter, Happy new 2020!
The steady state solve can be interpreted as zeroing the temperature field variation in time, ie, dT/dt = 0. The equation we use assumes velocity is 0 and is expressed with diffusivity, K, on the L.H.S, ie.
grad( Kgrad(T) ) = H/(rhocp)
So with H = 0 variations in rho don't affect the steady state equation because the R.H.S. = 0. Only variations of material diffusion will lead to a change in gradient of the geotherm.
In the model you sent via email the thermal gradients are discontinuous because the diffusivity values are also discontinuous, mantle.diffusivity != uppercrust.diffusivity.
Hope this helps.
— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or unsubscribe.
Sent from my iPhone
On 15 Jan 2020, at 9:20 am, Julian Giordani notifications@github.com wrote:
Hi Peter, Happy new 2020!
The steady state solve can be interpreted as zeroing the temperature field variation in time, ie, dT/dt = 0. The equation we use assumes velocity is 0 and is expressed with diffusivity, K, on the L.H.S, ie.
grad( Kgrad(T) ) = H/(rhocp)
So with H = 0 variations in rho don't affect the steady state equation because the R.H.S. = 0. Only variations of material diffusion will lead to a change in gradient of the geotherm.
In the model you sent via email the thermal gradients are discontinuous because the diffusivity values are also discontinuous, mantle.diffusivity != uppercrust.diffusivity.
Hope this helps.
— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub, or unsubscribe.
Hi Julian,
Happy New Year to you, too.
OK, I think that I understand the equation you're solving and why the solution is as it is. I'm still a bit confused how this solution achieves conservation of energy, but perhaps we should hold off for an in-person discussion with Patrice and Romain in February.
Roberta, Thanks for the plotting tips.
Cheers, Peter
Hi All,
I don't want to beat a dead horse, but have been thinking on this further and wish to raise the discussion once more.
Given my concerns with the discontinuities in heat flux coming from the UW steady state solution, I analytically solved the steady state heat equation in a rod consisting of 2 materials and no internal heat production, beginning with the formulation from Julian, grad (k grad(T)) = H / (rho cp), where k is thermal diffusivity: In doing so I realized that in addition to boundary conditions at the ends and continuity of temperature across the material boundary, an additional constraint is needed to relate the temperature gradient across the material boundary. UW (and many or most other dynamic solutions of the heat equation) appears to solve the problem such that k(grad(T)) is continuous across the boundary (where k is diffusivity). However, I don't think there is any physical basis for this constraint. Conservation of energy dictates that heat flux should be constant across the boundary; thus the additional constraint across the material transition should require that Kappa(grad(T)) is continuous across the material transition, where Kappa is thermal conductivity.
Below is a figure illustrating the potential effect on lithospheric geotherms. For illustrative purposes, I assume a two layer lithosphere in which crust and mantle have equivalent thermal conductivity and heat capacity, but density differs (2780 vs 3300 kg/m3) and hence so does thermal diffusivity. With homogeneous conductivity and no radiogenic heating, the geothermal gradient should be linear (black line), but the UW solution is in green. In this case, with a thick crust (~50 km), the result is ~65 degrees (K) difference in temperature at the Moho.
My impression is that the UW solution is consistent with many numerical solutions of the heat equation. Is continuity of k grad(t) implicit in numerical solutions that use thermal diffusivity? Is this just an accepted uncertainty in the solution? If so I think it should be acknowledged.
I have not considered the extension of this analysis to 2D or 3D, extension to the dynamic solution, or how a solution that enforces continuity of Kappa(grad(T)) might be implemented numerically. The objective of my 1D analysis is only to understand what should be the correct solution and what are its implications in this simple case.
Cheers, Pete
Hi Peter We start with the 1D conduction-advection heat equation (in which K is the conductivity): [cid:clip_image001.png]
Applying this equation to an infinite rod with no internal heat source (A = 0) and no heat advection (uz = 0) then: [cid:clip_image002.png]
This means that the evolution of temperature along the rod follows a linear trend. To get to T(z) we integrate twice and we get successively: dt/dz = a then T(z) = a.z + b In this last equation b is the temperature at z=0 (let’s call this temperature To), and a is the temperature gradient dT/dz. According to the Fourier’s law: dT/dz = Q/K (with Q the heat flow and K the conductivity). Hence: a = Q/K . Therefore, we get to: T(z) = Q/K . z + To This equation is valid for all z as long as the rod has homogeneous conductivity K.
Now let's assume a rod made of two sections (z0 to z1 and z1 to z2) with different conductivity Ki and Kj. In this case we have: For z ≤ z1: T(z) = Q / Ki . z + To For z1 < z < z2 : T(z) = (Q / Ki . z1 + To) + Q / Kj . (z- z1) We note that the first term (Q / Ki . z1 + To) is the temperature at z = z1. We also note that, because there is no internal heating, the heat flow Q is a constant throughout the rod. This implies that: Q = Ki . (dT/dz)i = Kj . (dT/dz)j … and since we have assumed that Ki ≠ Kj, then - to keep the heat flow constant - the temperature gradients (dT/dz)i and (dT/dz)j must also be different. In other terms, the break in the geotherm slope is necessary to obey the conservation of energy.
Let me know if this help.
Cheers
Patrice
On 6 Mar 2020, at 2:54 am, plovely notifications@github.com<mailto:notifications@github.com> wrote:
Hi All,
I don't want to beat a dead horse, but have been thinking on this further and wish to raise the discussion once more.
Given my concerns with the discontinuities in heat flux coming from the UW steady state solution, I analytically solved the steady state heat equation in a rod consisting of 2 materials and no internal heat production, beginning with the formulation from Julian, grad (k grad(T)) = H / (rho cp), where k is thermal diffusivity: [image]https://user-images.githubusercontent.com/30320218/75995265-c7ae2800-5ec1-11ea-9e2b-837e11d55e8f.png In doing so I realized that in addition to boundary conditions at the ends and continuity of temperature across the material boundary, an additional constraint is needed to relate the temperature gradient across the material boundary. UW (and many or most other dynamic solutions of the heat equation) appears to solve the problem such that k(grad(T)) is continuous across the boundary (where k is diffusivity). However, I don't think there is any physical basis for this constraint. Conservation of energy dictates that heat flux should be constant across the boundary; thus the additional constraint across the material transition should require that Kappa(grad(T)) is continuous across the material transition, where Kappa is thermal conductivity.
Below is a figure illustrating the potential effect on lithospheric geotherms. For illustrative purposes, I assume a two layer lithosphere in which crust and mantle have equivalent thermal conductivity and heat capacity, but density differs (2780 vs 3300 kg/m3) and hence so does thermal diffusivity. With homogeneous conductivity and no radiogenic heating, the geothermal gradient should be linear (black line), but the UW solution is in green. In this case, with a thick crust (~50 km), the result is ~65 degrees (K) difference in temperature at the Moho. [image]https://user-images.githubusercontent.com/30320218/75997261-b581b900-5ec4-11ea-893d-b9dd3a3f969d.png
My impression is that the UW solution is consistent with many numerical solutions of the heat equation. Is continuity of k grad(t) implicit in numerical solutions that use thermal diffusivity? Is this just an accepted uncertainty in the solution? If so I think it should be acknowledged.
I have not considered the extension of this analysis to 2D or 3D, extension to the dynamic solution, or how a solution that enforces continuity of Kappa(grad(T)) might be implemented numerically. The objective of my 1D analysis is only to understand what should be the correct solution and what are its implications in this simple case.
Cheers, Pete
— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://github.com/underworldcode/UWGeodynamics/issues/118?email_source=notifications&email_token=ADI3DZOT47DFICN5BZMAGPDRF7DLXA5CNFSM4HG62SB2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEN5ZSXA#issuecomment-595302748, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ADI3DZNZPLFO3QGNPAR36LLRF7DLXANCNFSM4HG62SBQ.
Thanks, Patrice. I think that your analysis is consistent with mine (though, the figures didn't come through).
However, in the model used to generate the steady-state geothermal gradient (second figure of my previous post), Ki = Kj, but the geothermal gradient from UW is kinked because diffusivities are different (ki≠kj). Thus my concern that steady-state heat flow is discontinuous in a model with no internal heat source.
Cheers, Peter
Hi Peter
I don’t think there is anything wrong with underworld...
The Fourier’s law is clear: Q = K dT/dz. Assuming no internal heating, then the heat flow Q is constant throughout the lithosphere.
In this case, a kink in the geotherm at the Moho is only possible when there is a difference in conductivity K (with K = diffusivity . Cp . density).
If the conductivity is constant, then there should be no kink no matter the respective values of diffusivity, Cp and density. In other terms: diffusivity, Cp and density can change but if K remains constant, then the geotherm is linear.
“Conservation of energy dictates that heat flux should be constant across the boundary; thus the additional constraint across the material transition should require that Kappa(grad(T)) is continuous across the material transition, where Kappa is thermal conductivity."
Kappa is not the conductivity, it is the diffusivity. Could this be the source of the confusion?
Cheers
Patrice
On 6 Mar 2020, at 10:43 am, plovely notifications@github.com<mailto:notifications@github.com> wrote:
Thanks, Patrice. I think that your analysis is consistent with mine (though, the figures didn't come through).
However, in the model used to generate the steady-state geothermal gradient (second figure of my previous post), Ki = Kj, but the geothermal gradient from UW is kinked because diffusivities are different (ki≠kj). Thus my concern that steady-state heat flow is discontinuous in a model with no internal heat source.
Cheers, Peter
— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://github.com/underworldcode/UWGeodynamics/issues/118?email_source=notifications&email_token=ADI3DZOHGTMBKOZQWCKSVQLRGA2KDA5CNFSM4HG62SB2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEN7DXPY#issuecomment-595475391, or unsubscribehttps://github.com/notifications/unsubscribe-auth/ADI3DZKRFC2BAUSZIWGQSP3RGA2KDANCNFSM4HG62SBQ.
Hi Patrice, I'll send my Jupiter notebook when I get to the office in the morning. Yes, the notation is confusing but I'm quite sure I have set up the UW model so that conductivity and heat capacity are uniform but density and diffusivity vary between layers. Cheers, Peter
On Thu, Mar 5, 2020, 6:34 PM Patrice Rey notifications@github.com wrote:
Hi Peter
I don’t think there is anything wrong with underworld...
The Fourier’s law is clear: Q = K dT/dz. Assuming no internal heating, then the heat flow Q is constant throughout the lithosphere.
In this case, a kink in the geotherm at the Moho is only possible when there is a difference in conductivity K (with K = diffusivity . Cp . density).
If the conductivity is constant, then there should be no kink no matter the respective values of diffusivity, Cp and density. In other terms: diffusivity, Cp and density can change but if K remains constant, then the geotherm is linear.
“Conservation of energy dictates that heat flux should be constant across the boundary; thus the additional constraint across the material transition should require that Kappa(grad(T)) is continuous across the material transition, where Kappa is thermal conductivity."
Kappa is not the conductivity, it is the diffusivity. Could this be the source of the confusion?
Cheers
Patrice
On 6 Mar 2020, at 10:43 am, plovely <notifications@github.com<mailto: notifications@github.com>> wrote:
Thanks, Patrice. I think that your analysis is consistent with mine (though, the figures didn't come through).
However, in the model used to generate the steady-state geothermal gradient (second figure of my previous post), Ki = Kj, but the geothermal gradient from UW is kinked because diffusivities are different (ki≠kj). Thus my concern that steady-state heat flow is discontinuous in a model with no internal heat source.
Cheers, Peter
— You are receiving this because you commented. Reply to this email directly, view it on GitHub< https://github.com/underworldcode/UWGeodynamics/issues/118?email_source=notifications&email_token=ADI3DZOHGTMBKOZQWCKSVQLRGA2KDA5CNFSM4HG62SB2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEN7DXPY#issuecomment-595475391>, or unsubscribe< https://github.com/notifications/unsubscribe-auth/ADI3DZKRFC2BAUSZIWGQSP3RGA2KDANCNFSM4HG62SBQ>.
— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/underworldcode/UWGeodynamics/issues/118?email_source=notifications&email_token=AHHKMWQFBTE6ERBHJPAX6NLRGBAI3A5CNFSM4HG62SB2YY3PNVWWK3TUL52HS4DFVREXG43VMVBW63LNMVXHJKTDN5WW2ZLOORPWSZGOEN7NSCQ#issuecomment-595515658, or unsubscribe https://github.com/notifications/unsubscribe-auth/AHHKMWQQ45TNZ34Y4H3GNFTRGBAI3ANCNFSM4HG62SBQ .
Hi All, Just checking in if you've had any further discussion on this since our email exchange with Patrice March 6-7. I hope everyone is well amidst the COVID crisis. Thanks.
Hi Romain,
Another issue related to the same model I just sent but with appropriate diffusivity (Model.diffusivity = 7.72e-7 * u.metre**2 / u.second).
The model has no radiogenic heat production, so in steady-state the heat flux should be constant everywhere below the land-air interface.
In this model, diffusivity is homogeneous, as is heat capacity, so thermal conductivity should be proportional to density (k=alpha rho Cp).
Because density of the mantle is greater than density of the crust, conductivity of the mantle should be greater and, given q=-k * grad(T), thermal gradient should be less in the mantle than in the crust. Instead, what I see (figure below) is a constant gradient in the temperature field and relatively greater heat flux in the mantle than in the crust.
Note: I calculate heat flux in ParaView by: 1) calculating the gradient of the temperature field using the "Gradient of Unstructured DataSet" filter 2) calculate heat flux using the calculator: q=grad(T) k = grad(T) alpha rho Cp, where grad(T) and rho come from paraview and alpha and Cp are constants.
Thanks, Peter