Dynamic ice melting

[cbegeman]

Implements new Monin-Obukhov Similarity Theory (MOST) method for an ice-ocean interface (McPhee, 1981; McPhee et al., 1987) and solves for heat and salt fluxes due to melting using the three-equation parameterization (Asay-Davis et al., 2016)

• Activated by setting namelist parameter most_method = 'mcphee'
• Requires that namelist parameter constant_flux_layer = 'top'
• Fluxes at the top boundary of momentum, heat and salt can be set at the initial timestep to values given by namelist file. At subsequent timesteps they are solved for by surface_layer_fluxes_mod module
• New namelist parameter for the drag coefficient drag_coeff
• New namelist parameters for distance or index count from the boundary to define "far-field" properties for MOST method k_offset and z_offset
• Implements TEOS10 polynomial solutions for freezing temperature and derivative of freezing temperature with respect to salinity

[cbegeman]

@xylar I haven't yet determined how to treat momentum fluxes. I've implemented McPhee's (1981) depth-dependent equation here. The potential issue is that if the scaled depth zeta is too large, then momentum fluxes go to 0 so that essentially the modeled flow doesn't feel the boundary. The image below shows the scale factor tau as a function of zeta

[cbegeman]

Note: branched from PR#34

[cbegeman]

@xylar I've been cleaning up this branch and found that one of the commits is more closely related to https://github.com/xylar/palm_les_lanl/pull/4. Should I reopen and amend that branch or create a new branch?

[xylar]

@cbegeman, once a branch is merged, it should just be deleted. So it won't be possible to re-open #4 and add to it. You should just create a new branch from master and git cherry-pick the commit you want onto that branch.

[cbegeman]

@xylar Ready for your review. I changed the base to https://github.com/xylar/palm_les_lanl/pull/34 for ease of viewing only melting-associated commits. Note that commit https://github.com/xylar/palm_les_lanl/pull/36/commits/398383a7a702717e1ec60b153ca930aa758d95e0 is associated with https://github.com/xylar/palm_les_lanl/pull/39, but I've kept it here for ease of testing.

[xylar]

I've implemented McPhee's (1981) depth-dependent equation here. The potential issue is that if the scaled depth zeta is too large, then momentum fluxes go to 0 so that essentially the modeled flow doesn't feel the boundary. The image below shows the scale factor tau as a function of zeta.

I can't remember if we discussed this. I definitely see the issue. McPhee's parameterization makes sense when the shear between the ocean and the ice isn't cause by melt (e.g. when it's cause by ocean currents with some other origin or wind stresses on sea ice). But this is precisely why I'm worried about the applicability of the McPhee parameterization in our case, where I expect buoyancy rather than Coriolis to be dominant near the ice base and buoyancy, not friction/drag to be dominant deeper down (the Jenkins 2016 argument). But I think using McPhee as a starting point is okay. If zeta ends up being so large for some reason (e.g. stratification) that we're in the tails of the exponential even in the top layers in the model, I think we really would need to rethink if it's the right model. It's probably worth writing out zeta somehow to make sure it's not crazy.

[cbegeman]

I've implemented McPhee's (1981) depth-dependent equation here. The potential issue is that if the scaled depth zeta is too large, then momentum fluxes go to 0 so that essentially the modeled flow doesn't feel the boundary. The image below shows the scale factor tau as a function of zeta.

I can't remember if we discussed this. I definitely see the issue. McPhee's parameterization makes sense when the shear between the ocean and the ice isn't cause by melt (e.g. when it's cause by ocean currents with some other origin or wind stresses on sea ice). But this is precisely why I'm worried about the applicability of the McPhee parameterization in our case, where I expect buoyancy rather than Coriolis to be dominant near the ice base and buoyancy, not friction/drag to be dominant deeper down (the Jenkins 2016 argument). But I think using McPhee as a starting point is okay. If zeta ends up being so large for some reason (e.g. stratification) that we're in the tails of the exponential even in the top layers in the model, I think we really would need to rethink if it's the right model. It's probably worth writing out zeta somehow to make sure it's not crazy.

Sounds good. I have all the output (i.e., L_O, friction velocity) I need to recompute zeta from existing runs so we can see what kinds of values it takes under different conditions. Let's chat about it in our next meeting.

[xylar]

I just realized that this PR is pointing to top_constant_flux_layer as the destination branch. This should be changed to master.

[cbegeman]

@xylar You can now test this branch using https://github.com/xylar/palm_les_lanl_viz/blob/master/ice_shelf/vis_run_short_long.py and namelist files located in https://github.com/xylar/palm_les_lanl_viz/tree/master/ice_shelf/test_cases I included the run scripts so you can get a sense for the number of processors and run time. I haven't finished running the long case, so I'm not sure how long that will take (in run time). Let me know if you have any questions about this.

[cbegeman]

I've implemented McPhee's (1981) depth-dependent equation here. The potential issue is that if the scaled depth zeta is too large, then momentum fluxes go to 0 so that essentially the modeled flow doesn't feel the boundary. The image below shows the scale factor tau as a function of zeta.

I can't remember if we discussed this. I definitely see the issue. McPhee's parameterization makes sense when the shear between the ocean and the ice isn't cause by melt (e.g. when it's cause by ocean currents with some other origin or wind stresses on sea ice). But this is precisely why I'm worried about the applicability of the McPhee parameterization in our case, where I expect buoyancy rather than Coriolis to be dominant near the ice base and buoyancy, not friction/drag to be dominant deeper down (the Jenkins 2016 argument). But I think using McPhee as a starting point is okay. If zeta ends up being so large for some reason (e.g. stratification) that we're in the tails of the exponential even in the top layers in the model, I think we really would need to rethink if it's the right model. It's probably worth writing out zeta somehow to make sure it's not crazy.

Sounds good. I have all the output (i.e., L_O, friction velocity) I need to recompute zeta from existing runs so we can see what kinds of values it takes under different conditions. Let's chat about it in our next meeting.

For the strong thermal driving case, the zeta values are between -0.05 and -0.1. The McPhee stress scaling seems appropriate here in that the flow will still "feel" the effects of friction despite strong stratification (i.e., momentum fluxes are still fairly large relative to their boundary values). I'll keep an eye on this in future simulations, but for now we're looking good.

[xylar]

@cbegeman, great! I'll take a look either tomorrow or Wednesday, and we can talk about it more at your meeting on Wednesday.

[xylar]

For the strong thermal driving case, the zeta values are between -0.05 and -0.1.

Great, thanks for checking that! I agree, that wounds perfectly reasonable.

[cbegeman]

@xylar Maybe you can think through with me why freezing conditions (over some portion of the domain) would lead the simulation to crash. The "crash" in this case is some horizontal velocity component near the boundary equaling NaN. The problem could arise earlier than the velocity term being assigned, but all the terms in the melt rate parameterization are real. The crash also doesn't happen at the first timestep where freezing occurs but takes a few timesteps (15s model time) to occur.

[cbegeman]

@xylar Maybe you can think through with me why freezing conditions (over some portion of the domain) would lead the simulation to crash. The "crash" in this case is some horizontal velocity component near the boundary equaling NaN. The problem could arise earlier than the velocity term being assigned, but all the terms in the melt rate parameterization are real. The crash also doesn't happen at the first timestep where freezing occurs but takes a few timesteps (15s model time) to occur.

@xylar This issue is resolved by imposing both lower and upper limits on Obukhov length L_O which effectively set the range of stability parameter eta_star from ~0.25 to 1. The lower limit on L_O ensures that shear stresses near the interface don't get too large in high melting regions. This is equivalent to values of zeta that are too large as we discussed above. The upper limit on L_O ensures that when slow melting or freezing occurs, the stress shape function is equivalent to that under neutral stratification. See commit 11c676d

[xylar]

Great, I think that makes a lot of sense. The upper limit (eta_star = 1) is, without question, the right one. The lower limit is presumably somewhat arbitrary (thus eta_star ~ 0.25, rather than an exact number). It would be good to see if our simulations are sensitive to that lower limit but probably not until we know more about what typical simulations look like.

Maybe we need to start making a table of parameters and potential ranges, and start ranking them in terms of how relevant they are. We should also non-dimensionalize the problem and explore the space in terms of the nondimensional parameters. Typically, this reduces the number of parameters to explore (maybe something to discuss next week).

Very exciting that you figured out a reason for the issue and a sensible solution!

[cbegeman]

@xylar Thanks for reviewing! I found one notational error in https://github.com/xylar/palm_les_lanl/pull/36/commits/7ee508e610e8d9c2ed10237212f3767da0b8e6fe which does not affect computation. The rest of the commits are unchanged. I'm ready to merge this.

[xylar]

It looks like you rebased after the merge of #42, which is good because the redundant commit doesn't seem to be here anymore. Go ahead and merge!