Update master with the developments that went into the neXtSIM v2 paper by einola · Pull Request #662 · nansencenter/nextsim

einola · 2025-01-23T09:33:40Z

This PR contains the most recent developments that made it into the neXtSIM v2 paper. This includes

Fixes to the penetrating shortwave scheme
Checks on reading in MLD
Updates for reading in forcing files
The melt pond scheme
The impact of atmospheric stability on ice-ocean drag

This addresses issue #634, by adding I (the penetrating short-wave radiation) as an output to IABulkFluxes and using this in the Winton and Semtner models.

The value of I_0 is quite uncertain it seems, but let's use the Winton (2000) value for now (30%).

I put the arguments in the wrong place in the IABulkFlux function, compared to the call. This is fixed. I also added t1 and t2 as output variables for the moorings.

Issue634 pen sw

The initial code written by me and Dani. It compiles, but it needs an option to turn on the meltpond scheme and mooring outputs.

It had a misspelling in its name - and more importantly, the permeability was incorrectly calculated.

It's called thermo.use_meltponds

We can now output "meltpond_volume", "meltpond_lid_volume", and "meltpond_fraction". All are averaged over the entire grid cell.

There were lots of minor problems with the initial implementation. Most of them are now fixed. The model runs and produces not unreasonable looking meltponds. The main remaining problems are - For some reason the meltpond volume recovers almost immediately after flushing - I get very deep meltponds in the oldest and most ridged ice. These can be deeper than the ice is thick and they never go away, because they don't flush and the lid becomes more than a meter thick and grows very slowly.

The isPermiable function now takes all temperatures into account and selects the Assur or Notz formulation, depending on temperature.

Some code clean-up of isPermeable. No functionality change or bug fixing.

A less convoluted and more obvious way of growing, melting, or forming a lid on top of a melt pond.

There's a risk that the ponds become unrealistically deep. The can even become deeper than the ice is thick. This patch prevents this, and (arbitrarily) limits the depth to 1/3 of the ice thickness by flushing excess water via the runoff variable. The limit can be thought out more carefully.

I was missing the melting term for the lid. This is now fixed.

In this version, the meltponds don't affect the freshwater balance (EMP). I also added a code to remove meltponds if the lid gets too thick (more then 30 cms).

I simply added a constant (tunable) albedo for the fraction of ponds without a lid (or a lid thinner than 5 cm - arbitrary number). This has a substantial effect in initial tests, but much more testing is needed.

This adds an alb_scheme == 4 option. It simply assumes constant snow, ice, and meltpond albedos. Something like alb_ice = 0.54, alb_sn = 0.83, and alb_pnd = 0.3 makes sense.

…T in alb scheme 4 + fixes

This adds a check on the MLD in checkFieldsFast, which throws an std::runtime_error if MLD is outside [0, 10 km]. The lower bound is quite strict (could also be a few mm) and the upper bound is just there because the code requires an upper bound. In reality we just want to check the lower bound.

…T in alb scheme 4 + fixes

We're still in the CICE spirit, but this includes some important fixes.

This way, we don't have to calculate the specific humidity too often. It also fits better within the code.

A correct calculation of the Obukov length this time! Some minor additional adjustments.

I multiplied the flux with cp + cpv. Not sure why, it should just be cp.

But the default is to do it. This is really only useful for debugging!

The new drag was not given correctly to the dynamics because I forgot to multiply with wind speed and atmospheric density(!)

This moves the calculation of some of the more complex constants out of the element loop to save time. Also adds the limiting length scale as an option, to put a sensible limit on the Obukov length (Richard says 100 m is good).

I missed a factor two when transcribing the equations. Fixing this actually improves the results!

typo

I managed to introduce a bug in the ocean drag calculation in a misguided attempt for minor optimisation. Remember Knuth!

I went through the constants yet again, and everything looks fine. Also changed the default Linv limit to something that makes more sense (and what I've used for my tests the whole time).

Mpnds4

It should be the Mellor and Kantha one, as this is implied in the paper we're writing.

Change the default value of melt_type

Stability

Change licence to MIT

tdcwilliams

some minor clarifications to add + typos

tdcwilliams · 2025-01-23T11:20:05Z

+                    longName = "Meltpond fraction";
+                    stdName  = "meltpond_fraction";
+                    Units    = "1";
+                    cell_methods = "area: mean";


is this area fraction, and is it fraction of the total grid area or just of the ice-covered area? Perhaps change the long name to "Meltpond area fraction"

ps if it is the fraction of ice-covered area use cell_methods = "area: mean where ice";

same for the other new melt-pond and drag variables

after reading the rest of the code the drags should be "area: mean where ice" but the melt pond variables can stay as "area: mean"

tdcwilliams · 2025-01-23T11:24:27Z

+            ("thermo.force_neutral_atmosphere", po::value<bool>()->default_value(false),
+             "Don't modify the neutral (default) atmospheric drag coefficient to take atmospheric stability into account.")
+            ("thermo.limiting_lengthscale", po::value<double>()->default_value( 1. ),
+             "A limit for the Obukov lenght (m).")


Suggested change

"A limit for the Obukov lenght (m).")

"A limit for the Obukov length (m).")

tdcwilliams · 2025-01-23T11:28:35Z

-            M_age_det[i] += dt*M_conc[i];
+            double w_age = old_conc <= 0? 0.: std::min(old_conc/M_conc[i], 1.);
+            // no ice => w_age = 0 => age = dt
+            // melting => old_conc > M_conc => w_age = 1


Suggested change

// melting => old_conc > M_conc => w_age = 1

// melting => old_conc > M_conc => w_age = 1 => age += dt

tdcwilliams · 2025-01-23T11:29:20Z

-            double del_vi_thick = M_thick[i] - old_vol;
+            w_age = old_vol <= 0 ? 0. : std::min(old_vol/M_thick[i], 1.);
+            // no ice => w_age = 0 => age = dt
+            // melting => old_vol > M_thick => w_age = 1


Suggested change

// melting => old_vol > M_thick => w_age = 1

// melting => old_vol > M_thick => w_age = 1 => => age += dt

tdcwilliams · 2025-01-23T11:31:57Z

+                psih = 2.*std::log(0.5*(1.+x*x));
+            }
+
+            // 4. The drag coefficients: ( \frac{k}{\ln{z/z_0} - \Pshi} )^2


Suggested change

// 4. The drag coefficients: ( \frac{k}{\ln{z/z_0} - \Pshi} )^2

// 4. The drag coefficients: ( \frac{k}{\ln{z/z_0} - \Psi} )^2

tdcwilliams · 2025-01-23T11:35:22Z

+    // double const roff = (0.85 - 0.7*M_conc[cpt]) ; Holland et al. 2012
+    M_pond_volume[cpt] += (1-roff)*availableWater*M_conc[cpt];
+
+    // Flush the pond if there's not enough ice. Skip everyting if there's no pond.


Suggested change

// Flush the pond if there's not enough ice. Skip everyting if there's no pond.

// Flush the pond if there's not enough ice. Skip everything if there's no pond.

tdcwilliams · 2025-01-23T11:38:54Z

+    D_pond_fraction[cpt] = std::sqrt(M_pond_volume[cpt]/dep2frac);
+
+    // Make sure there is no pond on snow
+    D_pond_fraction[cpt] = std::min(D_pond_fraction[cpt], 1.-hs/(hs+0.2));


if hs = 0, D_pond_fraction = 1?

No, but it can only happen that D_pond_fraction = 1 if hs = 0.

tdcwilliams · 2025-01-23T11:47:08Z

+    M_lid_volume[cpt] += delLidVolume;
+    M_pond_volume[cpt] -= delLidVolume;
+
+    // Remove lid if the pond is frozen solid or if it's too thick


should M_thick and/or M_h_young increase if pond freezes solid?

No, the pond volume is just an inactive tracer. It doesn't hold any mass.

einola · 2025-01-23T12:07:41Z

Because of branch protection (both master and develop are protected), it's best to merge this as it is and then fix everything in a separate PR.

Only changes to comments and metadata for output files.

Address comments from Tim on PR #662

Merge pull request #662 from nansencenter/develop

einola and others added 30 commits May 30, 2023 15:39

Fix issues with penetrating short-wave radiation into the ice

8b16ca5

This addresses issue #634, by adding I (the penetrating short-wave radiation) as an output to IABulkFluxes and using this in the Winton and Semtner models.

Change the default I_0 to Winton's value

5cf14b5

The value of I_0 is quite uncertain it seems, but let's use the Winton (2000) value for now (30%).

Merge remote-tracking branch 'origin/master' into issue634_pen_sw

9e6bcbf

Bug fix to previous commits

df9539b

I put the arguments in the wrong place in the IABulkFlux function, compared to the call. This is fixed. I also added t1 and t2 as output variables for the moorings.

Merge pull request #641 from nansencenter/issue634_pen_sw

895d5e7

Issue634 pen sw

Initial meltpond implementation

9c4f7fa

The initial code written by me and Dani. It compiles, but it needs an option to turn on the meltpond scheme and mooring outputs.

Fix the permeability function

7ec61a3

It had a misspelling in its name - and more importantly, the permeability was incorrectly calculated.

Add a swithc to turn on meltpond scheme

ed0a979

It's called thermo.use_meltponds

GridOutput for meltpond variables

8888d9b

We can now output "meltpond_volume", "meltpond_lid_volume", and "meltpond_fraction". All are averaged over the entire grid cell.

Improved permiability funciton

c3661a7

The isPermiable function now takes all temperatures into account and selects the Assur or Notz formulation, depending on temperature.

Cleaning up isPermeable

dea95fa

Some code clean-up of isPermeable. No functionality change or bug fixing.

Better handling of the lid formation

04679dc

A less convoluted and more obvious way of growing, melting, or forming a lid on top of a melt pond.

Fix surface energy balance

8cdfcd9

I was missing the melting term for the lid. This is now fixed.

Meltponds don't affect EMP

9f4d561

In this version, the meltponds don't affect the freshwater balance (EMP). I also added a code to remove meltponds if the lid gets too thick (more then 30 cms).

Add pond fraction to CCM3 albedo scheme

e576de3

I simply added a constant (tunable) albedo for the fraction of ponds without a lid (or a lid thinner than 5 cm - arbitrary number). This has a substantial effect in initial tests, but much more testing is needed.

update env_compile_intel_linux.bash after fram update

7fabfc4

should not change l_aerobulk

342ee8c

Adding modules + moved nextsim_intel

58017f9

Add a new albedo scheme (4)

ede8045

This adds an alb_scheme == 4 option. It simply assumes constant snow, ice, and meltpond albedos. Something like alb_ice = 0.54, alb_sn = 0.83, and alb_pnd = 0.3 makes sense.

Merge branch 'develop' into metponds2

29395bb

limit mpf to 0.5 of non-ridged ice area + snow albedo decreases with …

949d75e

…T in alb scheme 4 + fixes

adding a namelist for melt pond max fraction coverage on level ice

fce07dc

update env_compile_intel_linux.bash after fram update

5c13fdb

bugfix dynamics-type no_motion

a1627c7

removing changes in env_intel_linux from this branch

f0b0443

update env_compile_intel_linux.bash after fram update

15b4472

limit mpf to 0.5 of non-ridged ice area + snow albedo decreases with …

50a22d6

…T in alb scheme 4 + fixes

einola and others added 18 commits August 1, 2024 13:48

Fixes to the drag coeff scheme

eae5771

We're still in the CICE spirit, but this includes some important fixes.

Calculate drag coeffs in AIBulkFormula

e5d099f

This way, we don't have to calculate the specific humidity too often. It also fits better within the code.

Improvements to stability calculation

0ff5e7e

A correct calculation of the Obukov length this time! Some minor additional adjustments.

Bug fix for sensible heat flux

c4019f1

I multiplied the flux with cp + cpv. Not sure why, it should just be cp.

Make stability calculations optional

137d7a1

But the default is to do it. This is really only useful for debugging!

Fix drag in dynamics

574dc5a

The new drag was not given correctly to the dynamics because I forgot to multiply with wind speed and atmospheric density(!)

Cleaning of the Obukov code

c6e36ef

This moves the calculation of some of the more complex constants out of the element loop to save time. Also adds the limiting length scale as an option, to put a sensible limit on the Obukov length (Richard says 100 m is good).

Bug correction for psim stable

fe05a67

I missed a factor two when transcribing the equations. Fixing this actually improves the results!

Update gridoutput.hpp

14ba377

typo

Bug fix for ocean drag

f1c8677

I managed to introduce a bug in the ocean drag calculation in a misguided attempt for minor optimisation. Remember Knuth!

Final polishing of the stability alculators

79319c2

I went through the constants yet again, and everything looks fine. Also changed the default Linv limit to something that makes more sense (and what I've used for my tests the whole time).

Merge remote-tracking branch 'origin/develop' into mpnds4

15351ae

Changing location of namelist reading + cleaning

cafcd6b

Merge pull request #655 from nansencenter/mpnds4

2276129

Mpnds4

Merge remote-tracking branch 'origin/develop' into stability

716d74a

Change the default value of melt_type

f0ce868

It should be the Mellor and Kantha one, as this is implied in the paper we're writing.

Merge pull request #660 from nansencenter/melt_type

fd9348f

Change the default value of melt_type

Merge pull request #656 from nansencenter/stability

1882328

Stability

einola requested a review from tdcwilliams January 23, 2025 09:33

einola marked this pull request as draft January 23, 2025 09:34

Merge pull request #663 from nansencenter/master

9bcf49d

Change licence to MIT

einola marked this pull request as ready for review January 23, 2025 10:02

tdcwilliams requested changes Jan 23, 2025

View reviewed changes

einola requested a review from tdcwilliams January 23, 2025 12:08

tdcwilliams approved these changes Jan 23, 2025

View reviewed changes

tdcwilliams merged commit 9383d90 into master Jan 23, 2025

einola added a commit that referenced this pull request Jan 24, 2025

Address comments from Tim on PR #662

f1efe28

Only changes to comments and metadata for output files.

tdcwilliams added a commit that referenced this pull request Jan 26, 2025

Merge pull request #664 from nansencenter/pr662_minor_fixes

61c4f2b

Address comments from Tim on PR #662

einola added a commit that referenced this pull request Jan 28, 2025

Merge pull request #665 from nansencenter/master

7463e22

Merge pull request #662 from nansencenter/develop

	"A limit for the Obukov lenght (m).")
	"A limit for the Obukov length (m).")

	// melting => old_conc > M_conc => w_age = 1
	// melting => old_conc > M_conc => w_age = 1 => age += dt

	// melting => old_vol > M_thick => w_age = 1
	// melting => old_vol > M_thick => w_age = 1 => => age += dt

	// 4. The drag coefficients: ( \frac{k}{\ln{z/z_0} - \Pshi} )^2
	// 4. The drag coefficients: ( \frac{k}{\ln{z/z_0} - \Psi} )^2

	// Flush the pond if there's not enough ice. Skip everyting if there's no pond.
	// Flush the pond if there's not enough ice. Skip everything if there's no pond.

Conversation

einola commented Jan 23, 2025

Uh oh!

tdcwilliams left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

einola commented Jan 23, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

einola commented Jan 23, 2025 •

edited

Loading