Excessive velocity when including thermal model

clairewcook

I have been running a stress balance + mass transport model for an accumulating ice sheet, without the thermal model. I set it to have a constant temperature throughout the domain and define the rheology (B value) based on that constant temperature throughout the domain. For the rheology, I use parameters for grain boundary sliding accommodated basal slip (from Goldsby and Kohlstedt 2001). I have run this model with a temperature of 200 K throughout the domain (I'm modeling ice flow on Mars) and also with a temperature of 210 K throughout the domain. The maximum surface velocities reached in those models are ₀.4 m/yr and 0.8 m/yr respectively. However, I then tried to incorporate the thermal model. I set a constant surface temperature of 200 K and a basal heat flux of 30 mW/m². I modified one of the rheology law functions (nye.m) in order to use the rheology parameters I chose rather than the built in ones, and I set md.material.rheology_law='NyeH2O' which (I assume) uses that nye rheology law function. There are various issues that have come up when running this, including excessively high (and low) temperatures in some parts of the domain at late time steps (which may to be related to instabilities that I need to make other changes to solve). However, even before these unphysical changes in temperature, when the temperature throughout most of the ice is between 200 K and ₂₁₀ K, the maximum velocity in the ice reaches an unrealistic 40 m/yr, higher than the max velocity reached when the temperature is 210 K throughout the ice. I'd appreciate any ideas on why that may be.

mathieumorlighem

Hello @clairewcook
With the thermal model "on", you have to be very careful with boundary conditions. You will need Dirichlet conditions along all inflow boundaries. I recommend to set spctemperture not only at the surface but along all lateral boundaries to avoid advecting unconstrained temperature (e.g., pos = find(md.mesh.vertexonboundary); md.thermal.spctemperature(pos) = XX);).
Now for the NyeH2O rheology, I recommend taking a look at the C++ function itself: src/c/shared/Elements/NyeH2O.cpp. Is it consistent with what you had in mind?
I hope this helps
Mathieu

clairewcook

mathieumorlighem Thanks for your response. I had not thought to set Dirichlet conditions for temperature on all of the inflow boundaries. Just to clarify though, when spctemperature is not set for some boundary, is there is no default constraint defined on that boundary? And is it possible to set other types of boundary conditions besides Dirichlet, for example, if I wanted to set a no flux boundary condition along the inflow boundaries?

Regarding the rheology, I had just been using the MATLAB pre-compiled binaries for ISSM and modifying the bin/nye.m file, so I had not looked at the C++ function itself. It's good to know that's the important one to look at. I'm assuming now that if I wanted to define a custom rheology, I would need to get the source code, change that C++ function, and then recompile ISSM with that change?

Thanks,
Claire

mathieumorlighem

Hi @clairewcook

If you don't set spctemperature, there is no default value for the temperature and your problem may be ill-posed at the inflow boundary (i.e. where the ice velocity is going inward).

You can impose zero flux (∇T.n = 0), you actually don't have to do anything because that will be a right hand side = 0, BUT you need to be careful because this will only work for the diffusion term. The advecting term (v.∇T) may not be 0 if the velocity is not 0 and so then you may be advecting in an unknown amount of heat. So, in short, if you want to impose zero flux, do not set spctemperature along the boundary but you need to make sure the speed is 0 (you will need to spcvx and vy and vz to 0). I hope this makes sense!

and yes, if you want to define a new custom rheology, you will need to change the C++ code and recompile. It's not complicated at all.
Don't hesitate if you need more help!
Mathieu