0-layer sea-ice model in SCM #46
Conversation
|
|
||
| # Semtner's Zero-layer Model | ||
|
|
||
| Sea ice is approximated as a 0-layer slab. This means that the model does not occupy a particular gridpoint between the mixed layer (ocean) and the atmosphere. It is essentially a column-wise ODE, with no horizontal transport between the columns. In the absence of ice, the sea-ice model reduces to the slab ocean formulation. The ice is assumed to have a negligible heat capacity, so there is no heat storage within it. The only store of energy comes from the ice thickness. |
There was a problem hiding this comment.
This is a little unclear to me - "The ice is assumed to have a negligible heat capacity, so there is no heat storage within it. The only store of energy comes from the ice thickness." What is the difference between heat storage and energy storage?
There was a problem hiding this comment.
I guess this is just saying that the ice layer doesn't have a prognostic equation for T. I will rephrase.
| $$ | ||
|
|
||
| ## 3. Transition from no ice to ice | ||
| - When the updated $T_{ml}^{t+1} < T_{melt}$, set $T_{ml}^{t+1} = T_{melt}$ and grow ice via $h_i$ due to the corresponding energy deficit. |
There was a problem hiding this comment.
Does this mean that when ice is present, T_ml quickly becomes T_freeze (due to F_base) and then remains there (F_base = 0)?
There was a problem hiding this comment.
It does :) (unless ocean_qflux is switched on)
| $$ | ||
| - while ice-covered conditions require that: | ||
| $$ | ||
| \rho_w c_w h_{ml}\frac{dT_{ml}}{dt} = - F_{base} |
There was a problem hiding this comment.
This is the same sign (-F_base) as in the ice equation. they should be opposite signs, right?
There was a problem hiding this comment.
Sorry, which ice equation do you mean? All fluxes should be positive when pointing upward, following the original paper.
There was a problem hiding this comment.
Ahh, you mean h_ice? The sign should be the same - e.g. an upward (positive) flux will decrease ocean temperature (T_ml) as well as decreasing the ice thickness (h_ice). Lmk if this is not what you meant. 🤔
| - When the updated $T_{ml}^{t+1} < T_{melt}$, set $T_{ml}^{t+1} = T_{melt}$ and grow ice via $h_i$ due to the corresponding energy deficit. | ||
| - When the updated $h_i^{t+1} <= 0$ from a non-zero $h_i^t$, adjust $h_i^{t+1} = 0$ and use the surplus energy to warm the mixed layer. | ||
|
|
||
| ## 4. Surface temperature ($T_s$) |
There was a problem hiding this comment.
If there is ice, T_s = surface temp of ice? Does that contradict the statement earlier that the ice is an isothermal slab? It seems like we have ice with surf temp = T_s, bottom temp = T_freeze, and flux at surface = F_i = k_i \frac{T_{melt} - T_s}{h_i}, flux at bottom = F_base, but we assume F_atm = F_i.
When the surface is ice free, is T_s = T_ml? Should mention that here?
There was a problem hiding this comment.
Yeah, I the phrasing wasn't quite right, due to an earlier plan that didn't materialize. Good catch! ;)
| dT/dt = ∇ μ ∇ T | ||
| where | ||
| T = 280 K at z = zmax_atm | ||
| dT/dt = - ∇ F_sfc at z = zmin_atm |
There was a problem hiding this comment.
Is this a lower BC? I'm confused that it seems to be setting dT/dt. Do we want -\mu \nabla T = F_sfc at the bottom?
|
|
||
| ## set BCs | ||
| bcs_bottom = Operators.SetValue(Geometry.WVector(F_sfc ./ p.p.λ)) | ||
| bcs_top = Operators.SetValue(Geometry.WVector(FT(0))) |
There was a problem hiding this comment.
Does not match BC in comments (T = 280)?
| dT/dt = - ∇ F_sfc at z = zmin_atm | ||
|
|
||
| We also use this model to calculate and accumulate the downward surface fluxes, F_sfc: | ||
| F_sfc = - λ * (T_sfc - T1) |
There was a problem hiding this comment.
Im a little confused why F_sfc is not F_rad+SHF + LHF here.
There was a problem hiding this comment.
ah OK, I see that this is a to-do item, and that F_sfc is consistent with the ice model (but not as described in the docs for the slab ice model). Maybe you could note in the documentation that as a stand-in for F_sfc, you use this expression?
| end | ||
|
|
||
| # update state | ||
| @. Y.T_ml += ΔT_ml |
There was a problem hiding this comment.
Does this work? I thought you could only change prognostic variables via callbacks (or at least, that it was only recommended to do so).
Would it be possibly to implement the ice model as a diff eq now, rather than implement this kind of funky "set Y by timestepping by hand, wrap inside a ODE problem which doesnt ever change dY from 0" (also, is dY always 0 when initialized? perhaps you should still set dY .= zero field inside the ice tendencies stub)
There was a problem hiding this comment.
Yeah, it works lol. I did it this way mainly to deal with the time stepping and to be able to access the ice state from previous time steps for diagnostics. We can do something more sophisticated in the next iteration, but it may not be straightforward since you need to do stuff between time integrating h_ice and T_sfc.
| @. ΔT_base = p.T_freeze - T_sfc | ||
| end | ||
| # surface is ice-covered, so update T_sfc as ice surface temperature | ||
| @. Y.T_sfc += ΔT_base |
There was a problem hiding this comment.
another item to keep in mind (we need to do so on the land side too) is that ice/ocean and snow/soil have different albedos and latent heats involved (when computing LHF). They may also have different roughness lengths.
| ΔT_base = similar(F_atm) | ||
|
|
||
| # ice thickness and mixed layer temperature changes due to atmosphereic and ocean fluxes | ||
| if h_ice[1] <= 0 # ice-free |
There was a problem hiding this comment.
Im not sure if this really matters, but Im not sure how often h_ice will be =0. You could say h_ice <eps(FT) instead to account for tiny difference from zero.
Do you see h_ice as negative ever?
There was a problem hiding this comment.
It shouldn't be negative and the @. Δh_ice = -h_ice should take care of this, but since it is not an integer I left it there in case there were some FP errors. I can invert the argument and change it to h_ice[1] > 0 for ice covered and else otherwise?
| T_top = FT(280.0), # fixed temperature at the top of the domain_atm [K] | ||
| ## slab ice parameters | ||
| ρc_ml = FT(4e6), # for mixed layer | ||
| F_basal = FT(120), # W / m2 / K |
There was a problem hiding this comment.
this is minor, but F_basal for the coefficient and F_base for the flux is a little confusing. maybe even a letter besides "F" for the coefficient would be helpful, since "F" denotes a flux everywhere else?
There was a problem hiding this comment.
Yeah, I agree. Great point!
| end | ||
|
|
||
| # adjust if transition to ice-covered | ||
| if (T_ml[1] + ΔT_ml[1] < p.T_freeze) |
There was a problem hiding this comment.
If h_ice > 0, I assume this can still be satisfied. That might be what you want, in any case! To keep T_ml at T_freeze and creat more ice. But thought Id point it out.
There was a problem hiding this comment.
Umm, h_ice > 0 wouldn't indicate a transition, so it would loop over the entire ice-covered domain, rather than just where it's needed. Thanks for pointing it out though!
| # adjust if transition to ice-free | ||
| if ((h_ice[1] > 0) & (h_ice[1] + Δh_ice[1] <= 0)) | ||
| @. ΔT_ml = ΔT_ml - (h_ice + Δh_ice) * p.L_ice / (p.h_ml * p.ρc_ml) | ||
| @. Δh_ice = -h_ice |
There was a problem hiding this comment.
It doesn't need it here. Once we can deal with spatially varying BCs in ClimaCore, this will be point wise, so all the index brackets should go away. :)
| return | ||
| end | ||
|
|
||
| function ∑tendencies_ice_stub(du, u, p, t) |
There was a problem hiding this comment.
Are you required to use forward Euler with this model, to handle the ice <-> no ice transitions?
| dY = du | ||
| FT = eltype(dY) | ||
|
|
||
| solve_ice!((; u = u, p = p), p.Δt) # timestepping outside of DeffEq (which is only used for logging here) |
There was a problem hiding this comment.
If you cant get around this, is it preferable to not make the ice variables prognostic (needing dY...) and instead just make them aux vars that get updated?
There was a problem hiding this comment.
I'm still using DiffEq to save the solution with time. Not sure it's easily possible otherwise.
| @. ΔT_base = p.T_freeze - T_sfc | ||
| end | ||
| # surface is ice-covered, so update T_sfc as ice surface temperature | ||
| @. Y.T_sfc += ΔT_base |
There was a problem hiding this comment.
the += here could be an issue. What if in the previous step, Y.T_sfc = mixed layer temp? In this case, would this give you the right ice surface temp? It seems safer to use Y.T_sfc = ice surface temp.
There was a problem hiding this comment.
This method should account for this. The ice surface temp = Y.T_sfc + ΔT_base (btw the name ΔT_base was ambiguous so I changed it).
There was a problem hiding this comment.
Hi Lenka. This generally seems fine, though the documentation is ahead of where the code is (in terms of F_sfc, and in terms of writing the ice/ocean equations as continuous ODEs rather than discretized equations).
Right now, it's a hard-coded algorithm for updating T_sfc, T_ml, h_ice using foward Euler stepping. If it will stay this way, I think it might be better to remove these from Y and stick them in p, because the code isnt actually timestepping them with ODE.jl/making use of the prob and keywords a user might set...
a couple items flagged that seemed like they could lead to issues - modifying the state directly, using a += for the ice covered T_sfc, and not setting dY = 0 each step (since you are modifying them without using dY - if it is nonzero at initialization, would ODE algorithm update Y?)
|
Thanks so much for all the above comments, @kmdeck. The docs are a little ahead in view of the next stage of this issue (#42 ), but you're right, it is clearer to say that explicitly. I tried to note where we're using simplifications for the time being (e.g. |
hi @LenkaNovak! Yes, feel free to merge! It looks like it is working based on your checks. And I'm not sure how to better handle the timestepping either - we'll face the same issue with snow on land. It makes sense to do what works now, and we can address it later, as you say. |
sea ice + conservation bc fix conservation check incl. ice clean up Kat fixes format
|
bors r+ |
|
Build succeeded: |
Steps
Proposed assignees
Proposed Reviewers
TODO in following PRs