Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

For Antarctica simulations in CESM, SMB can accumulate in the ocean #39

Closed
billsacks opened this issue Aug 25, 2021 · 34 comments
Closed
Labels
bug Something isn't working

Comments

@billsacks
Copy link
Member

I noticed that, in some of my Antarctica tests in CESM, SMB accumulates in areas that are in the ocean. I have noticed this in an I compset test that uses very coarse land resolution (10 deg x 15 deg: /glade/scratch/sacks/ERS_Vnuopc_Ly3.f10_f10_ais8_mg37.I1850Clm50SpGa.cheyenne_intel.20210811_223534_2a9bfb). It's possible that this problem won't arise in more realistic production resolutions (it doesn't appear in a T compset where the land forcing is at f09_g17 resolution; I haven't tried an I compset at that 1-degree resolution), but I still think that this is a real issue that should be fixed.

Here is a figure showing SMB from CISM's history file from the above test:

20210825_154540_6UE7AI

This positive SMB in the periphery of the grid (which I think should be open ocean) leads to an accumulation of ice thickness in these grid cells.

Here is the corresponding SMB map from the coupler history file:

20210825_163521_wvZ5Rz

From digging in to some of the coupling fields and from comparing the behavior with that in Greenland simulations, I think what's going on is:

  • The icemask correctly starts as 0 in open ocean points in Antarctica.
  • There is always potentially a non-zero SMB passed in areas outside of the ice sheet: CTSM generates SMB outside as well as inside the icemask, and assumes that CISM will discard any SMB sent outside the icemask.
  • For Greenland runs, CISM appears to properly be discarding SMB sent outside the icemask. However, for Antarctica runs, there are some grid cells outside the icemask that accept this SMB and start growing ice.
  • It appears that these problematic points start with usurf = 0. So it seems like the problem is that some usurf = 0 points are able to accept SMB in Antarctica, in contrast to Greenland. Interestingly, there is a swath of points in the Antarctica run that seems to zero out the SMB, but beyond that swath, it accepts SMB.

This is a problem not only because it leads to ice growth in the open ocean, but also because I think it would break conservation in a fully-coupled run: CTSM assumes that the icemask (actually icemask_coupled_fluxes, but in practice they are effectively the same) dictates where CISM is able to receive SMB. Conversely, it assumes that, if icemask_coupled_fluxes is 0, then CISM is not able to receive SMB there, and so CTSM should send any generated SMB directly to the ocean, via the snow capping flux. But, for system-wide conservation, that means that CISM needs to be consistent in this: If CISM says that icemask_coupled_fluxes is 0 (which is typically the case over open ocean), then it needs to discard any SMB that is sent there. I'll admit that this coupling is rather subtle and error-prone, and could probably use to be reworked (probably with the coupler/mediator handling some of this logic itself), but for now that's how things work.

@billsacks billsacks added the bug Something isn't working label Aug 25, 2021
@gunterl
Copy link
Contributor

gunterl commented Aug 26, 2021 via email

@billsacks
Copy link
Member Author

@gunterl - interesting. I was attributing the blockiness to the coarse resolution land forcing (10° x 15°), but I hadn't looked into it carefully. Your explanation could explain why there is a ring of correctly-handled points in the ocean grid cells near Antarctica and the problems appear further afield.

The case is here: /glade/scratch/sacks/ERS_Vnuopc_Ly3.f10_f10_ais8_mg37.I1850Clm50SpGa.cheyenne_intel.20210811_223534_2a9bfb

@whlipscomb
Copy link
Contributor

@billsacks, Thanks for the detailed explanation and plots. Like @gunterl, I noticed that the SMB patterns seem connected to CISM's inactive blocks – the regions that are identified by CISM's ice_domain_mask as permanently ice-free, so that no processors need to be assigned there. We typically use this setting for Antarctica but not for Greenland.

My interpretation of the problem is that you're seeing a nonzero SMB in the inactive blocks. This might be an initialization issue. The coupler is passing a nonzero SMB over the ocean, and CISM fails to zero it out because it's not doing any local computations in those regions. Maybe we need to do something with a global array. These regions do exist at the global level, even though they aren't mapped to any local processor.

The inactive-block option is enabled by setting compute_blocks = 1 in the namelist. So a first step would be to set compute_blocks = 0 and see if this SMB issue goes away.

Then we can look for the best place in the code to add some new logic. I don't know where this would be, but I'm glad to take a look.

@billsacks
Copy link
Member Author

It looks like compute_blocks is already 0 in /glade/scratch/sacks/ERS_Vnuopc_Ly3.f10_f10_ais8_mg37.I1850Clm50SpGa.cheyenne_intel.20210811_223534_2a9bfb/run/cism.ais.config.

This is definitely something we should look at carefully when doing a run that has compute_blocks = 1, but that doesn't seem to be the issue here.

@whlipscomb
Copy link
Contributor

@billsacks –  Interesting. I also see you're running Antarctica on only 72 processors, which isn't enough to generate this pattern of active and inactive blocks.

Nevertheless, I think this problem is related to the compute_blocks option, because the pattern is so similar to what I'm used to seeing in Antarctic 8km runs with this option turned on. Could you please point me to your CISM input file? Maybe this file has a bad 'topg' field from a previous run with compute_blocks = 1.

@whlipscomb
Copy link
Contributor

@billsacks – Thinking ahead a bit, we'll want to make sure that when we do set compute_blocks = 1, we write out topg with values that won't mislead the coupler on restart. I think we want to avoid topg = 0 in the inactive regions. In these regions, we may want to save the initial global topg field and always write out this field, instead of a field that is gathered from the active blocks and has incorrect values for the inactive blocks.

@billsacks
Copy link
Member Author

The init file is: /glade/p/cesmdata/cseg/inputdata/glc/cism/Antarctica/ISMIP6_Antarctica_8km.1950_init_c210613.nc

@billsacks
Copy link
Member Author

billsacks commented Aug 26, 2021

And you're right: here's the topg from that file (the brownish points are exactly 0):

image

@billsacks
Copy link
Member Author

Besides the obvious issue that topg should have reasonable values for the whole domain instead of having 0 values in these more distant ocean points, this makes me think that there may be a more subtle issue of how points with usrf exactly equal to 0 are treated.

For the sake of determining what points are inside the icemask, we use this function:

logical function is_in_active_grid(geometry, i, j)
! Return true if the given point is inside the "active grid". The active grid includes
! any point that can receive a positive surface mass balance, which includes any
! point classified as land or ice sheet.
type(glide_geometry), intent(in) :: geometry
integer, intent(in) :: i, j ! point of interest
real(dp) :: usrf ! surface elevation (m)
! TODO(wjs, 2015-03-18) Could the logic here be replaced by the use of some existing
! mask? For now I am simply re-implementing the logic that was in glint.
usrf = thk0 * geometry%usrf(i,j)
if (usrf > 0.d0) then
! points not at sea level are assumed to be land or ice sheet
is_in_active_grid = .true.
else
is_in_active_grid = .false.
end if
end function is_in_active_grid

The assumption there (if I remember correctly) is that points with usrf <= 0 will throw away any SMB they are handed. But here we have points with usrf exactly equal to 0 that seem to be accumulating SMB. Is it possible that there is a check somewhere else in CISM that throws away SMB for points with usrf strictly less than 0, and that we should bring those two checks into sync so that points with usrf exactly 0 are treated consistently?

@whlipscomb
Copy link
Contributor

@Bill Sacks – Sorry, I should have looked at this file more closely. It's not an ideal input file. It seems to be a restart file that was repurposed as an input file, which is why topg is wrong. It also contains some fields that are needed only for exact restart and not for initialization.

In the near term, I think it will be enough to fix the topography. @gunter Leguy, could you please make a file that has the correct bedmachine topg throughout the domain?

In the past, it wasn't easy to remove files from inputdata subdirectories. Is this still true? It would be nice to be able to experiment with inputs before we settle on a production version.

I'll follow up with a reply to your most recent post.

@gunterl
Copy link
Contributor

gunterl commented Aug 26, 2021 via email

@billsacks
Copy link
Member Author

In the past, it wasn't easy to remove files from inputdata subdirectories. Is this still true? It would be nice to be able to experiment with inputs before we settle on a production version.

Once you have added a file to the inputdata svn repo, it is effectively frozen. But you can experiment with files on glade to come up with something you're happy with, then we can just add the final version.

Ideally the new input file will use exactly the same grid as the old one – then we won't need to remake any grid / mesh files.

@whlipscomb
Copy link
Contributor

@billsacks – Yes, CISM may be doing the wrong thing where topg is exactly zero. For instance, here in glad_timestep:

      ! Set acab to zero for ocean cells (bed below sea level, no ice present)                                                                            
      where (GLIDE_IS_OCEAN(instance%model%geometry%thkmask))
         instance%acab = 0.d0
      endwhere

In general, I've been deprecating the old Glide masks in favor of Glissade masks, so I'm not happy to see GLIDE_IS_OCEAN here. However, we would get the same result using the Glissade ocean mask. Over in glide_masks.F90, we have this:

!Identify points where the ice is floating or where there is open ocean                                                                                  
where (topg - eus < (-rhoi/rhoo) * thck)
    mask = ior(mask, GLIDE_MASK_OCEAN)   ! GLIDE_MASK_OCEAN = 8                                                                                          
elsewhere
    mask = ior(mask, GLIDE_MASK_LAND)    ! GLIDE_MASK_LAND = 4                                                                                           
endwhere

Cells with topg = thck = 0.d0 will be identified as land, hence capable of receiving an SMB.

I'd like to (1) bring these two checks into agreement, and I'd prefer to keep the convention that topg = 0 implies land, not ocean. In addition, we need to (2) fix the topg output for compute_blocks = 1, and (3) produce a better input file. Maybe you could do (1), I could work on (2), and in due course @gunterl could take care of (3).

@whlipscomb
Copy link
Contributor

Going forward, let's experiment with input files on Glade. These will be on the same 8km grid as the file that's already in inputdata. We won't put any more files in inputdata without some thorough vetting.

I anticipate that a 4km grid will be the default for production runs.

@whlipscomb
Copy link
Contributor

@gunterl – For now, I think it will be enough to give @billsacks a file with the correct topography and ice thickness. The thickness is already OK, since we have thk = 0 in the inactive region. So just swap out the topography.

But we should start thinking about fields that should go in the final input file. We can start with the fields that would be in a standard ISMIP6 file before spin-up. To these we would add tempstag, powerlaw_c_inversion, and deltaT_basin from the end of the spin-up, since both these fields evolve during the spin-up. I can't think of any more that we need, assuming that we have in mind a forward run with spun-up temperature, C_p, and deltaT_basin.

To prepare the final version, we should repeat our spin-up with compute_blocks = 0, so we don't have the block structure anywhere in the input file.

@billsacks
Copy link
Member Author

This sounds like a good plan. For (1), though, I'm not clear on exactly what you envision to keep these checks in agreement. The easiest thing would be for me to change the conditional in is_in_active_grid to check usrf >= 0 instead of usrf > 0. Is that what you have in mind for now, or do you feel that I should (a) use an existing mask variable in place of that check (if so, which one), and/or (b) change the GLIDE_IS_OCEAN check in glad_timestep to instead use a mask in glissade (if so, again, which one?).

@whlipscomb
Copy link
Contributor

@billsacks – What I had in mind is just that you change the conditional in is_in_active_grid to check usrf >= 0. I'll take care of GLIDE_IS_OCEAN at a later stage, probably as part of an overall mask cleanup. It's OK for now since it agrees with the glissade mask ('ocean_mask') that would replace it.

Thinking some more about compute_blocks = 1: This is a useful feature for people who want to do multiple long spin-ups in TG mode, but the cost savings is in the noise for century-scale BG runs. For now, it will be easier to work with input files that are defined on the entire grid. Down the road, to support compute_blocks = 1, would it make sense to write special values in the inactive-block region to avoid confusion?

Another issue: It doesn't feel quite right to me to put input files with spun-up temperature and basal friction fields permanently in inputdata. If we tweak any inversion parameter, we get a different answer. There's an analogous issue with spun-up BGC fields in other components, including CTSM. Are multiple files with spun-up results put in inputdata, or just the basic information used for a cold start? For ice sheets, the basic data include ice thickness and topography, geothermal heat flux, and climatological thermal forcing. These data could serve as the starting point for many spin-ups and wouldn't need to be replaced so often.

@gunterl
Copy link
Contributor

gunterl commented Aug 27, 2021 via email

@billsacks
Copy link
Member Author

@whlipscomb - Okay, thanks, I'll make that simple change shortly. Regarding your other questions:

Down the road, to support compute_blocks = 1, would it make sense to write special values in the inactive-block region to avoid confusion?

That makes sense to me.

Another issue: It doesn't feel quite right to me to put input files with spun-up temperature and basal friction fields permanently in inputdata. If we tweak any inversion parameter, we get a different answer. There's an analogous issue with spun-up BGC fields in other components, including CTSM. Are multiple files with spun-up results put in inputdata, or just the basic information used for a cold start? For ice sheets, the basic data include ice thickness and topography, geothermal heat flux, and climatological thermal forcing. These data could serve as the starting point for many spin-ups and wouldn't need to be replaced so often.

For CTSM, we put fully spun-up files in the inputdata repository, because this is what gets pointed to out-of-the-box so that people can run with a spun-up land and not need to do the expensive and somewhat tricky spinup themselves. These spun-up files are roughly 1 GB in size, with some of them a few GB (see /glade/p/cesmdata/cseg/inputdata/lnd/clm2/initdata_map). Even for a given version of the model, there are a number of different initial conditions files available (probably at least 10), for different years, atmospheric forcings, etc. My sense is that they are replaced every few years, when there is a new scientifically-supported version of the model. A few 10s of GBs is actually relatively small in the grand scheme of the inputdata repository: For example, generating all new surface datasets for CTSM generates a few TB of new data.

Following that example, I feel it would be very reasonable for you to include your current best estimate of inverted fields in the initial file, if you feel that's what users should use for their own simulations in general. But if this is something that you change frequently (say, multiple times a year), then it could make sense to update the out-of-the-box file – which would be the file in inputdata – less frequently, e.g., timed with releases of CISM.

@gunterl - thank you for adding the file. However, the format needs to be changed: Some of our supported systems have problems reading NetCDF4-formatted files, so we aren't even allowed to add NetCDF4-formatted files to the inputdata repository. I don't have experience with the built-in NetCDF compression: does this require NetCDF4 format? If not, you could simply convert this to a different format (in the past I have used ncks -5 FILENAME to convert to an allowed format). However, unless you know for sure that compressed files can be read by CESM on all systems, can you please send an email to Jim Edwards (cc'ing me) describing how you are doing the compression and confirming that it will be okay to do this? (It may be relevant to add to him that CISM doesn't use PIO.)

billsacks added a commit to billsacks/cism that referenced this issue Aug 31, 2021
As per suggestions from Bill Lipscomb in
ESCOMP#39, I am treating usrf == 0 as
land rather than ocean. This makes the conditional here consistent with
similar checks elsewhere in CISM.
billsacks added a commit that referenced this issue Aug 31, 2021
Make usrf condition in is_in_active_grid consistent with other checks

As per suggestions from Bill Lipscomb in
#39, I am treating usrf == 0 as
land rather than ocean. This makes the conditional here consistent with
similar checks elsewhere in CISM.
@billsacks
Copy link
Member Author

As I discussed with @whlipscomb by phone this afternoon, it didn't work to change the conditional in is_in_active_grid to just be usrf >= 0 instead of usrf > 0: All ocean points have usrf == 0, so this change made the icemask equal to 1 everywhere, which is not what we want.

@whlipscomb had some thoughts about how to handle this, such as checking usrf > 0 .or. topg >= 0, but since we want to be careful to get this right (and consistent with other CISM code) and since this is somewhat tangential to the main Antarctica issue here, we are going to defer this. I'll open a separate issue for this problem.

billsacks added a commit to billsacks/cism that referenced this issue Sep 1, 2021
…r checks"

This reverts commit 3b899f4.

This didn't work: it resulted in the ice mask being 1 everywhere.

See ESCOMP#39 (comment)

We'll fix this later (see ESCOMP#41)
@billsacks billsacks mentioned this issue Sep 1, 2021
billsacks added a commit that referenced this issue Sep 1, 2021
Revert PR #40

This didn't work: it resulted in the ice mask being 1 everywhere.

See #39 (comment)

We'll fix this later (see #41)
@whlipscomb
Copy link
Contributor

@billsacks and @gunterl – I just sent a detailed email. Here, I'd like to summarize my understanding of where we are and ask some questions.

  • Bill S. is now working with an input file that has several unnecessary fields, some of which aren't defined over the whole grid because they came from a run with compute_blocks = 1. At this stage, we want to run with compute_blocks = 0.
  • We can either (1) do a quick patch to that file, putting in a correct topg field, and hope that the missing regions in other fields don't mess things up, or (2) step back and make a new input file that more resembles what we would run with later.
  • If (2), then we can choose between (2a) a standard input file used for cold start, or (2b) an input file with a spun-up ice thickness, temperature and other fields.
  • It may be easier to prepare (2a) than (2b). But I've forgotten why we chose a spun-up run to start with. Bill S., for the work you're doing now, do you know whether it matters?

@billsacks
Copy link
Member Author

Thanks for your thoughts on this @whlipscomb . For the work I'm doing now, all I really need is something that can run for a few years and give results that aren't complete garbage. So probably any of (1), (2a) or (2b) will work. If you feel like it's worth the time to go ahead with (2a) or (2b) now, I can wait for those. But if you have other higher priority things you'd like to be doing now, I think it would also work to use (1) for now and then return to this later: I know it means putting a temporary file in inputdata that we plan to replace, but this file isn't too large.

@whlipscomb
Copy link
Contributor

@billsacks, I suggest (2a). Earlier this week, I was thinking (1) because I didn't want to hold you up. But now I'm thinking that we've never run the model with compute_blocks = 0 after reading input fields from a run with compute_blocks = 1 (nor would we, in practice). This seems a bit dangerous – we might spend time tracking down errors that result from an artificial configuration. I can work with @gunterl on a robust cold-start input file.

@whlipscomb
Copy link
Contributor

@billsacks, To get you going quickly (I hope), please look at this directory: /glade/p/cesm/liwg_dev/lipscomb/cism_antarctica/stage_cism_spinup

I'm using an input file that @gunterl created for ISMIP6 Antarctic spin-ups:
ISMIP6_Antarctica_8km.bedmachine2.nc

This config file is similar to what we used for ISMIP6, with a few updates to account for recent code changes:
ismip6.config.init

To make sure things are working, I launched a 1000-year run using the cism_driver executable compiled from cism_main_2.01.004:
1e376e1 (tag: cism_main_2.01.004, main) Merge pull request #36 from ESCOMP/lipscomb/ais_coupled

And the output looks OK. I think this will give you what you need for testing, after tweaking the namelists. Please let me know if you have any questions.

@billsacks
Copy link
Member Author

@whlipscomb and @gunterl - This new file (ISMIP6_Antarctica_8km.bedmachine2.nc) has lat and lon values that differ by more than just roundoff from the file we had been using before (the files in /glade/p/cesmdata/cseg/inputdata/glc/cism/Antarctica). Can one of you help reconcile this? Ideally the new file would use the same grid as the original one to avoid needing to update the mesh file pointed to by cime, but if we need to update the mesh file then we can do that. I can't tell, though, whether the differences in lat & lon values are intentional or if (for example) this new file is meant to be on the exact same grid as the old ones and the lat and lon values are just incorrect on this new one.

Specifically, I'm seeing:

  • The longitude convention differs between the new and old files (0 to 360 vs. -180 to 180); because of that I haven't looked closely for any possible differences in longitude values.
  • Latitude values look similar, but seem to differ by more than roundoff. For example, the latitude values of the first few grid cells are -57.8708610534668, -57.9235382080078, -57.9761734008789 on the old file, and -57.82823, -57.88089, -57.93351 on the new file.
  • We want latitude and longitude values to be double precision; this is the case on the old file but not the new file.

Could one of you please look into this and let me know how to proceed? The ideal solution for me would be if you determine that the new file is supposed to have the same latitude and longitude values as the old one, and then we can just copy the latitude and longitude values from the old one to the new one. But I realize that may not be the correct thing to do.

@whlipscomb
Copy link
Contributor

@billsacks, I will defer these questions to @gunterl. What you're calling the 'new' file was actually made some time ago for ISMIP6 simulations, and what you're calling the 'old' file may be a file for which Gunter computed lat and lon more recently. My guess would be that the 'old' values are the ones we want, but Gunter can confirm.

@billsacks
Copy link
Member Author

@whlipscomb - As I started to compare the config file you pointed me to with what comes out of the box, I remembered that this is tricky because there are quite a few intentional differences between what we use in CISM-in-CESM vs. these standalone configurations, which I spent a while going through a few months ago. So instead I dug up the config file you pointed me to a few months ago, which I used as the basis for the defaults in CISM-wrapper, and diffed this new config file with that one.

I am attaching the two files for reference.

ismip6.config.init_20210903.txt
ismip6.config.init_20210602.txt

Potentially relevant differences are the following:

  • temp_init: was 4, your new config file has 3. I see a note in the namelist xml file: "temp_init = 3 can only generate an advective-diffusive-balance profile if acab is present in the input file (which is currently only the case for the 4km input file). But it defaults to a linear profile (temp_init 2) if acab is not present. So to keep things simple, we specify temp_init = 3 for all input files, rather than making this dependent on glc_grid." It looks like acab is not on the Antarctica init file you pointed me to, so (if that comment is still correct) I guess that means it will use the fall-back behavior of a linear profile. Just confirming that's what you want here.
  • adjust_input_thickness: was false, your new config file has true
  • which_ho_cp_inversion: was 2, your new config file has 1
  • which_ho_bmlt_basin_inversion: was 2, your new config file has 1
  • Your new config file has a bunch of inversion_* parameters that weren't present in the old config file you gave me. These are all available in the namelist xml file, but with default values that differ from what you sent. I can update the defaults to what you sent for Antarctica, but wanted to confirm that these are relevant and therefore worth updating.

Can you please confirm that I should make all of the above changes?

@whlipscomb
Copy link
Contributor

@billsacks – Thanks for looking carefully at the new config file. Let me respond to each point:

  • temp_init = 3: The namelist comment is not quite accurate. We can get an advective-diffusive profile with temp_init = 3 if either acab or smb is present in the input file. (The only difference between acab and smb is the units.) So we should modify the comment.
  • In the new config file, I set adjust_input_thickness = true to patch up an issue with Lake Vostok in the input topg and thk fields.
  • For spin-up, we don't have a pre-existing powerlaw_c_inversion field, so we need to invert for this field during spin-up. Hence which_ho_cp_inversion = 1. Likewise for bmlt_basin inversion.
  • Yes, let's update the inversion defaults consistent with the values in the new config file.

In short: Yes, please make all the above changes. Thanks!

@billsacks
Copy link
Member Author

Thanks @whlipscomb - I'll make those changes soon. One more question about the inversion_* parameters: Is it safe to change the defaults for both gris & ais (which would be slightly easier), given that we use which_ho_cp_inversion and which_ho_bmlt_basin_inversion = 0 by default for gris? Or should we maintain the old default values for these inversion_* parameters for Greenland?

@whlipscomb
Copy link
Contributor

@billsacks – Yes, it's OK to change the defaults for both ice sheets. I haven't run Greenland with bmlt_basin inversion, but for Cp inversion the same parameters that work for the AIS should work for the GrIS.

@gunterl
Copy link
Contributor

gunterl commented Sep 9, 2021 via email

billsacks added a commit to ESCOMP/CISM-wrapper that referenced this issue Sep 9, 2021
See some comments in ESCOMP/CISM#39 for
details. (The main intent of the new init file is to address that issue,
but there are other reasons for updating the init file as well.)
@billsacks
Copy link
Member Author

Thanks for your recent input, @whlipscomb and @gunterl . I have started testing with the new init file and config settings.

@billsacks
Copy link
Member Author

Changes are here: ESCOMP/CISM-wrapper#62

@billsacks
Copy link
Member Author

With the new Antarctica init file, the problem that started this issue is fixed. I will make a CISM-wrapper tag with this change shortly.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working
Projects
None yet
Development

No branches or pull requests

3 participants