Fix routing/runoff for coupled model applications

[sanAkel]

Statement:

Routing table file (e.g., CF0180x6C_TM1440xTM1080-Pfafstetter.TRN) that gets generated by make_bcs script(s) puts run off in wrong locations:

• In the middle of ocean (tile) points (some could be islands?).
• On land (tile) points.

See following figure, locations are marked in panels (c) and (d). Attached is a terrain map of the Indonesian seas (W Pac) region that can be compared to panel (d). These have major consequences for the coupled model: causing it to crash and leaks in the hydrological budget, to name a few. Comments from @atrayano, @rdkoster and @gmao-rreichle will help.

@yvikhlya has been running this utility to fix above problem to some extent.

History:

• @amolod has asked for a fix since many years.
• More recently effort to centralize boundary conditions (files): @biljanaorescanin and generation (source code) by @weiyuan-jiang et al., requires us to take a fresh look at it.

List of steps to generate a routing file that is currently used:

For e.g., /discover/nobackup/projects/gmao/ssd/aogcm/atmosphere_bcs/Icarus-NLv3/MOM6/CF0180x6C_TM1440xTM1080_newtopo/CF0180x6C_TM1440xTM1080-Pfafstetter.TRN

1. Run make_bcs package. This will produce a TRN file, but it directs runoff to the land, that is incorrect (see below figures).

2. Take the TRN file from above step 1 and run @yvikhlya's reroute.py utility (see code at https://github.com/GEOS-ESM/GMAO_Shared/blob/main/GEOS_Util/coupled_diagnostics/g5lib/src/reroute_runoff/reroute.py) which fixes the routing so it ends up being in the ocean.

3. Incorporate @lcandre2 's runoff pathways so that runoff from glaciers is routed to the ocean in a more realistic way. Where is the script/code? ETA: I have this code (though is now several years old). It assigns tiles to a 'topographically appropriate' ocean cell. Because the current routing scheme simply puts runoff in the geographically nearest ocean cell (which is incorrect as frequently as it is correct). The longer term fix - which needs to be started - is to develop catchments and a pfafstetter routing scheme for landice and treat landice as a true child of surface.

4. Possibly (if the coupled model or ocean component blows up at run time) iterate the above 2 and/or 3 if too much freshwater is sent to the ocean. This will be done manually by modifying the ocean mask used in steps 2 and 3, such modified ocean mask, used for routing should be saved somewhere.

Note:

• Step 4 cannot be codified.
• Current make_bcs script(s) does not include above steps 2- 4, they are being done outside (manually) as described ⬆️ .

Action items:

• We need to have steps 2 and 3 included in the make_bcs script(s).

@yvikhlya 's comments on how to get :arrow_up: done:

Actually, the problem is in this utility

https://github.com/GEOS-ESM/GEOSgcm_GridComp/blob/develop/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/mk_runofftbl.F90

What it should do at the very start is to read tile file, read MOM mask and mark all tiles which fall on MOM land as land (change code from 0 to 20). 
Then do what it currently does. This would eliminate re-routing step.

If manual manipulation with routing is necessary, we would feed this utility a fake mask. Done.

cc: @wmputman, @tclune, @atrayano, @weiyuan-jiang

[sanAkel]

Thank you @gmao-rreichle.

I'm eternally hopeful! :)

[zyj8881357]

Some updates on my side about this issue:

The higher-res outlet files I sent to Atanas show the final outlet locations (for each cell) in lon and lat indexes (not degrees) on a global 1-arcmin raster grid. All the outlet locations are on land, so they must be linked to ocean points for the use in the coupled model.

As I checked the data and mask information, the problem is that there may some "lakes" (water bodies) laying between the land and ocean. The runoff should go through these "lakes" to ocean. If some of these "lakes" are not explicitly simulated in the ocean model, it will be not easy to link the outlet on land to a correct ocean point (for example, simply linking the nearest ocean point to the outlet on land).

But the good thing is that, for all the ~22000 ocean outlets in the world, most of them only deliver a small amount of runoff to ocean. In these cases, a small shift of the outlet location in ocean will not be a big problem I guess? So we may use the nearest ocean point as the outlet in ocean for these cases. For all other outlets (mean discharge > 100 m3/s, considering that the first 150 major rivers in the world have a mean discharge >= 2000 m3/s), I can calculate the distances between the outlet on land and the nearest ocean point. If the distance is very small (<300 m), it is also OK for linking them to the nearest ocean point? If the distance is not small, I can do manually to link this outlet to an ocean point. As I checked, only ~60 outlets needs to be set manually. Not a big deal.

So if I have a land-ocean mask used in the ocean model, I can do all this above. But as I learnt from the former posts, the land-ocean mask can be changed when we change the ocean model (MOM, MIT). So a better way may be that, on my side, I use some high-res mask file (eg. the 10s GEOS land-lake-ocean mask file) to determine the lat and lon degrees (in ocean) for each outlets with the methods above. And then others can link this degree information to ocean tiles very easily (for example, using the nearest ocean tile) no matter what land-ocean mask is in use.

[yvikhlya]

Each ocean grid has its own land/sea mask and there are no lakes on any of ocean grids. Ocean mask (MOM/MIT) is not the same as GEOS or realistic mask. Runoff should be routed from land tile into ocean tile which corresponds to wet grid cell in ocean model.

[atrayano]

@zyj8881357 Thanks for your long update. The way I see it, we need to produce ocean mask file(s), and typically those vary by the ocean model and the ocean resolution. I like your idea to "locale" the nearest ocean outlet point using a high-res definition. This will be an improved version of the current scheme. But you do this process only once. Then we run mk_runofftbl.x (which I am supposed to modify and read/process the relevant ocean mask)

[zyj8881357]

@zyj8881357 Thanks for your long update. The way I see it, we need to produce ocean mask file(s), and typically those vary by the ocean model and the ocean resolution. I like your idea to "locale" the nearest ocean outlet point using a high-res definition. This will be an improved version of the current scheme. But you do this process only once. Then we run mk_runofftbl.x (which I am supposed to modify and read/process the relevant ocean mask)

Yes. In a short words, I can provide you files with outlet locations in lat and lon degrees (rather than in indexes). And this version should be very easy to link to the ocean tiles in ocean models.

[atrayano]

Let me check the assumptions in the current version of the code. I am nearly 100% sure that it expects that lat-lon info (either in degrees or as index) points to ocean cell (not land cell)

[zyj8881357]

Let me check the assumptions in the current version of the code. I am nearly 100% sure that it expects that lat-lon info (either in degrees or as index) points to ocean cell (not land cell)

Yes, in the current version (that I mentioned above, not the version I already sent to you), the lat-lon info points to ocean (not land). The ocean is defined in the mask file of "GEOS5_10s_mask.nc". There may still some difference for the ocean mask between "GEOS5_10s_mask.nc" and ocean model (MOM/MIT), but the (ocean->ocean shift) will be much easier to handle than before (land->ocean shift). If you think this way is OK for you, I will go to proceed this updated version and send it to you.

[sanAkel]

@zyj8881357 thank you!

Can you upload an ascii file with:

River (name)	lat (deg N)	lon (deg E)
ABCD	y_1	x_1

You should be able to upload it here, if not, please share via google drive; assuming no access to discover cluster at NCCS.

Thanks again.

[zyj8881357]

@zyj8881357 thank you!

Can you upload an ascii file with: River (name) lat (deg N) lon (deg E) ABCD y_1 x_1

You should be able to upload it here, if not, please share via google drive; assuming no access to discover cluster at NCCS.

Thanks again.

I uploaded it to the google drive. The link is: https://docs.google.com/spreadsheets/d/1FbQsqoYL-knSmF_qzKlOcC8Pp3kxW4Ok/edit?usp=share_link&ouid=110154394974887000797&rtpof=true&sd=true

It's a spreadsheet. I also added the mean discharge for each river. Please check whether these locations are consistent with the ocean mask in the ocean model.

[sanAkel]

Thanks @zyj8881357 I downloaded the file. (You can remove it from above URL if you wish to.)

Question: Is river(catchID) unique? In other words, does it change with tile file resolution?

[zyj8881357]

Thanks @zyj8881357 I downloaded the file. (You can remove it from above URL if you wish to.)

Question: Is river(catchID) unique? In other words, does it change with tile file resolution?

Thanks for the downloading! The catchment IDs are unique and do not change with the tile file resolution. In other words, the outlets locations in the table do not change with the tile resolution.

[sanAkel]

Thanks for the downloading! The catchment IDs are unique and do not change with the tile file resolution. In other words, the outlets locations in the table do not change with the tile resolution.

In that case,

1. Say the resolution of the atmosphere is c12, which is very coarse, how do you propose to deal with more than one catchment ID having the same (lon, lat)?

2. See attached figure/notebook, what happened to Greenland? You don't include glaciers?

   

[gmao-rreichle]

@sanAkel:

Say the resolution of the atmosphere is c12, which is very coarse, how do you propose to deal with more than one catchment ID having the same (lon, lat)?

I presume this would have to be handled in the mk_runofftbl.F90, and it's not different from what would have been needed with the older outlet info. What Yujin provides is the best location of the outlet based on the currently used Pfafstetter basin/watershed classification, which has a fixed "resolution".

[...] what happened to Greenland? You don't include glaciers?

Yujin only provides the outlets from the land. I don't know how exactly @lcandre2 handles glacier runoff from Greenland. There should be some (non-glaciated) land along the coast, especially in the north of Greenland. In your map above, I don't see outlets from Greenlands "land". That's something for Yujin to look into. Possibly Yujin and Lauren left all of Greenland to Lauren's processing?

[zyj8881357]

Thanks for the downloading! The catchment IDs are unique and do not change with the tile file resolution. In other words, the outlets locations in the table do not change with the tile resolution.

In that case,

1. Say the resolution of the atmosphere is `c12`, which is _very_ coarse, how do you propose to deal with more than one catchment ID having the same (lon, lat)?

2. See attached figure/notebook, what happened to Greenland? You don't include glaciers?

Thanks for the feedback! Now I added three outlets for the Greenland (as shown by the attached figure). These three outlets are for the land area (non-glaciated) in the Greenland (as they have river IDs: 154779, 154795 and 161496). I did not treat the glacier area because we don't have the river network information (ie, Pfaf code) associated with the glacier area. You can access the updated version of the outlet locations at: https://docs.google.com/spreadsheets/d/1rOG4KhQLbP4YJKYhqP6fqIkHoG3wjz5w/edit?usp=share_link&ouid=110154394974887000797&rtpof=true&sd=true

As mentioned by @gmao-rreichle, this table only lists the locations of all outlets to ocean. It cannot be directly used in the code (ie, read_riveroutlet.F90 and mk_runofftbl.F90) for making bcs. At this stage, please just check if all this points are not difficult for you to transfer to the ocean area in your ocean models. If some of them (especially for those with mean discharge > 100 m3/s) are not easy for you to transfer to your ocean area, just tell me their river IDs and show me the land-ocean mask you are using. So I can fix it further. After all this outlet locations are OK with you to use, I will go to proceed them to the format that can be easily read in the read_riveroutlet.F90 and mk_runofftbl.F90. Thank you!

[gmao-rreichle]

Thanks for the update, @zyj8881357. Re. Greenland, I wasn't sure if we have Pfafstetter info for the non-glaciated land area. It sounds like we do not. Are the three outlets you added for Greenland meant to collect runoff from all of the non-glaciated land in Greenland? This doesn't seem right to me. Shouldn't there be a much larger number of outlets surrounding the entire coast? I suppose this should be looked at in conjunction with how @lcandre2 routes the glacier runoff towards the coast.

[lcandre2]

Thanks for the tag @gmao-rreichle. @zyj8881357, we do not currently have pfastetter codes for Greenland or Antarctica because most runoff flows through the ice sheet and then along the ice sheet bed to a land or ocean outlet. I have a routing code that addresses this and provides better output locations for greenland and antarctica (in the rare cases it has it) runoff. However, it is offline. @zyj8881357, can we touch base via a teams meeting to coordinate how to approach this? It might be the right time to develop catchments for the ice sheets. We have the data, I just need to do the work.

[zyj8881357]

Thanks for the update, @zyj8881357. Re. Greenland, I wasn't sure if we have Pfafstetter info for the non-glaciated land area. It sounds like we do not. Are the three outlets you added for Greenland meant to collect runoff from all of the non-glaciated land in Greenland? This doesn't seem right to me. Shouldn't there be a much larger number of outlets surrounding the entire coast? I suppose this should be looked at in conjunction with how @lcandre2 routes the glacier runoff towards the coast.

If you look at the figure, the green area is for the glacier area and there is no Pfafstetter code associated to it. In our GEOS5_10arcsec_mask.nc file, there are only three unit catchments around it. Their Pfafstetter code tells us they are isolated catchments without any upstream and downstream catchments. I agree with you that our Pfafstetter code is too simple over the Greenland, and it would be better for us to treat the Greenland in conjunction with glacier people.

[lcandre2]

Thanks for the update, @zyj8881357. Re. Greenland, I wasn't sure if we have Pfafstetter info for the non-glaciated land area. It sounds like we do not. Are the three outlets you added for Greenland meant to collect runoff from all of the non-glaciated land in Greenland? This doesn't seem right to me. Shouldn't there be a much larger number of outlets surrounding the entire coast? I suppose this should be looked at in conjunction with how @lcandre2 routes the glacier runoff towards the coast.

If you look at the figure, the green area is for the glacier area and there is no Pfafstetter code associated to it. In our GEOS5_10arcsec_mask.nc file, there are only three unit catchments around it. Their Pfafstetter code tells us they are isolated catchments without any upstream and downstream catchments. I agree with you that our Pfafstetter code is too simple over the Greenland, and it would be better for us to treat the Greenland in conjunction with glacier people.

@zyj8881357 We must have been writing at the same time! see my message above. Let's chat.

[zyj8881357]

Thanks for the tag @gmao-rreichle. @zyj8881357, we do not currently have pfastetter codes for Greenland or Antarctica because most runoff flows through the ice sheet and then along the ice sheet bed to a land or ocean outlet. I have a routing code that addresses this and provides better output locations for greenland and antarctica (in the rare cases it has it) runoff. However, it is offline. @zyj8881357, can we touch base via a teams meeting to coordinate how to approach this? It might be the right time to develop catchments for the ice sheets. We have the data, I just need to do the work.

Yes. I just saw it! We can set up a teems meeting to see how to coordinate the land and glacier. If you think it is better to discuss the meeting schedule via Email, you email me at yujinz@umbc.edu. Thank you!

[sanAkel]

@zyj8881357, Can I ask you for a help?

Could you make a diff plot of ALL land data from here which has SRTM?

To be precise, plot of land data from GEOS5_10s_mask.nc - ⬆️ URL. If the land difference is zero or what @GEOS-ESM/land-team think are sensible, for once, we can switch to a consistent source of topography for entire earth! 🚀

(of course, we can try/test first!)

[zyj8881357]

@zyj8881357, Can I ask you for a help?

Could you make a diff plot of ALL land data from here which has SRTM?

To be precise, plot of land data from GEOS5_10s_mask.nc - ⬆️ URL. If the land difference is zero or what @GEOS-ESM/land-team think are sensible, for once, we can switch to a consistent source of topography for entire earth! 🚀

(of course, we can try/test first!)

Hi Santha, I think you want me to compare the elevation diff between our land source data and the DEM from the link? I'd like to do the comparison, but the GEOS5_10s_mask.nc only provides the land-sea mask and does not include the elevation values. So I may need to get the elevation source data in our land model first before doing the comparison. @gmao-rreichle @biljanaorescanin

[sanAkel]

Hi Santha, I think you want me to compare the elevation diff between our land source data and the DEM from the link? I'd like to do the comparison, but the GEOS5_10s_mask.nc only provides the land-sea mask and does not include the elevation values. So I may need to get the elevation source data in our land model first before doing the comparison. @gmao-rreichle @biljanaorescanin

@zyj8881357 I recall, sometime ago when I plotted GEOS5_10arcsec_mask_freshwater-lakes.nc https://github.com/GEOS-ESM/GEOSgcm_GridComp/blob/5938f19b7acfb46ea44ad42f421448ee03762241/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs#L437 it had elevation over land and none (=0) over the open ocean. Maybe I'm wrong (I do forget things)! But, yes, elevation diff is what I'm looking for ...

I'll also look at ${GLOBAL_CATCH_DATA} that I'm feeding to the make_bcs. Information such as this has made me curious about any differences!

Thanks!

[biljanaorescanin]

@zyj8881357 I will send you our topo file. You can see more about that in README file you got with boundary conditions. Any version of BCs files you have will have under clsm directory a README file with more description on what your boundary conditions are and where they come from.

[gmao-rreichle]

@sanAkel, @zyj8881357: FYI, I think the topo file for the atmospheric model is smoothed, so even if the topo data underpinning the Pfafstetter watershed delineations and the atm topo data are both based on SRTM, I wouldn't expect a perfect match. And then there's the fact that SRTM is limited to ~60N/S and high-lat topo is from some other source. @biljanaorescanin

[sanAkel]

@sanAkel, @zyj8881357: FYI, I think the topo file for the atmospheric model is smoothed, so even if the topo data underpinning the Pfafstetter watershed delineations and the atm topo data are both based on SRTM, I wouldn't expect a perfect match. And then there's the fact that SRTM is limited to ~60N/S and high-lat topo is from some other source. @biljanaorescanin

Thank you @gmao-rreichle, I would like to see the differences and later ask you all ( @GEOS-ESM/land-team ) for all the steps you have taken to generate the topography. I expect to be able to generate it (elevation) myself, once I'm given directions ... of course, will not clobber that with this issue!

[biljanaorescanin]

I could also pickup AGCM topo files but I will need specific resolution you want to compare, since data is organized by resolution. These AGCM topo file are created by Ben A. maybe he has a raster for it I am not sure how that part works. For example: NCAR_TOPO_GMTED/TOPO_CF0720x6C/smoothed/topo_DYNave*.data So this is CF720:

[sanAkel]

@biljanaorescanin how did you read that binary topo_DYN_ave_*.data? Can you copy/paste your script/code? Thanks!

[biljanaorescanin]

Santha each file has .bin version and .nc4

/discover/nobackup/projects/gmao/bcs_shared/legacy_bcs/Icarus-NLv3/NCAR_TOPO_GMTED/TOPO_CF0720x6C/smoothed/gmted_DYN_ave_720x4320.nc4 this is one we link into topo_DYN_ave_720x4320.data

[zyj8881357]

Hi Santha, I think you want me to compare the elevation diff between our land source data and the DEM from the link? I'd like to do the comparison, but the GEOS5_10s_mask.nc only provides the land-sea mask and does not include the elevation values. So I may need to get the elevation source data in our land model first before doing the comparison. @gmao-rreichle @biljanaorescanin

@zyj8881357 I recall, sometime ago when I plotted GEOS5_10arcsec_mask_freshwater-lakes.nc

https://github.com/GEOS-ESM/GEOSgcm_GridComp/blob/5938f19b7acfb46ea44ad42f421448ee03762241/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/make_bcs#L437

it had elevation over land and none (=0) over the open ocean. Maybe I'm wrong (I do forget things)! But, yes, elevation diff is what I'm looking for ...

I'll also look at ${GLOBAL_CATCH_DATA} that I'm feeding to the make_bcs. Information such as this has made me curious about any differences!

Thanks!

Hi Santha, here is the comparison between the elevations used in the land model (srtm30_withKMS_2.5x2.5min.nc) and from your link (SRTM15+). While there is no substantial difference between them (panels a and b), but small random differences occur outside of 60S-60N (panel c). @sanAkel @biljanaorescanin @gmao-rreichle

[sanAkel]

Thank you @zyj8881357 for this update.

Locations (such as the Canadian coast) in the Arctic, Greenland, Antarctic, Western Pacific, they all seem to have differences- in other words, uncertainty in coastline is evident.

I always knew of such problematic spots when comparing our land-model mask with that from the ocean model (derived from GEBCO). Your plots clearly show that we can trace those differences to the source:

• Land model (probably some version of SRTM?).
• Ocean model (GEBCO 2021).

[gmao-rreichle]

@sanAkel, @zyj8881357 : The GEOS catchments were defined (and mean elevations derived) using SRTM data where available, ie., ~60S-60N. Outside of these lats, SRTM did not provide data. (The shuttle didn't fly poleward of 60 N/S lat.) I don't know what the "link to SRTM15+" is, but I suspect that the data at this link were backfilled outside of 60 N/S differently from what was done in the Raster processing. I'm not sure what this is comparison is meant to accomplish. I don't see a cause for concern. cc: @rdkoster @biljanaorescanin

[lcandre2]

Thanks for the tag @gmao-rreichle. @zyj8881357, we do not currently have pfastetter codes for Greenland or Antarctica because most runoff flows through the ice sheet and then along the ice sheet bed to a land or ocean outlet. I have a routing code that addresses this and provides better output locations for greenland and antarctica (in the rare cases it has it) runoff. However, it is offline. @zyj8881357, can we touch base via a teams meeting to coordinate how to approach this? It might be the right time to develop catchments for the ice sheets. We have the data, I just need to do the work.

All, @zyj8881357, @biljanaorescanin, apologies for the long delay. I have produced outlets and associated drainage basins for Greenland that are gridded to the GEOS510arcsec_mask.nc. The steps closely, but not completely follow that of Verdin (2013, Final Report: High Resolution Topographic Analysis for GMAO’s Catchment LSM) and use ArcticDEM v3 100m mosaic cropped to Greenland (Porter et al., 2018; https://www.pgc.umn.edu/data/arcticdem/). 100 m is approximately 10 arc seconds at a latitude of 70° N. Steps 3-6 were completed using Topotoolbox v2 (Schwanghart & Scherler, 2014; https://topotoolbox.wordpress.com/).

1.     Crop DEM Greenland only. 2.     Fill a small number of missing grid cells using bilinear interpolation. No extrapolation was used, and this step was only designed to remove NaN values assigned to a handful of locations across Greenland. 4.     Hydrologically condition DEM by filling sinks. 5.     Calculate flow direction and flow accumulation grids using the D8 routing algorithm. 6.     Delineate ‘surface streams’ and drainage basins using an area threshold of 250 km2, the same threshold as in Verdin (2013). Streams are extracted to vector format and contain the Strahler stream order numbering system; drainage basins retain a raster format. See below for streams and 7.     Extract drainage basin outlets. For Greenland, at the given resolution there are 525 outlet locations. 8.     Map drainage basins and outlets to the GEOS10arcsec_mask.nc grid using ‘nearest’. This mapping follows the methodology outlined in Verdin (2013) & Mahanama et al. (2015; GMAO Tech Report 39). 9.     Confirm each point is appropriately assigned a drainage basin number.   Notes (in no particular order) 1.     @zyj8881357 expect an email about this and about how you would like the data transferred shortly. 2.     I have done this for all of Greenland (both ice and land) because Yujin indicated a lack of outlets across Greenland and because previous efforts excluded Greenland. I can also produce drainage basins and outlets for ice (or ice sheet) only points as well. This would reduce the number of outlet locations around Greenland by excluding many of the islands. 3.     Outlet locations are on land and still need to be moved to the nearest ocean cell. 4.     This information should be sufficient to create the routing table. While Pfaffstetter interbasins can be determined using the surface streams from step 5, including Greenland ice sheet Pfaffstetter coding system would likely cause issues (or require a rearrangement of the tile numbering) in subsequent processes because the 200000000 ice value would be replaced. My suggestion would be to append add a basin number to the Greenland ice values instead – though this is completely up for discussion because I do not know what the exact impacts would be. 5.     In the original routing scheme, is that water from land ice (200000000) locations is routed to the nearest ocean location using an expanding spiral search, so some modification of the code will be needed to exclude the Greenland ice points because my previous work modified the routing table after the spiral search. 6.     There should be future consideration of how to deal with glacier runoff routing. 7.     Because there is no dynamic ice sheet, water is theoretically routed across the surface. In reality, much of the water produced on the ice surface only flows short distances before flowing englacially to the ice sheet bed. The traditional routing scheme for subglacial water accounts for the overburden pressure of the ice and basal topography in the potential gradient calculation; however, to substantially differ from surface routing the bed slope must exceed the surface slope by ~10x, which is uncommon. 8.     Though the streams and drainage basins reach the ice sheet divide, runoff generally only occurs along the ice sheet edges.

Note: middle panel legend refers to stream order.

[sanAkel]

Thanks for the tag @gmao-rreichle. @zyj8881357, we do not currently have pfastetter codes for Greenland or Antarctica because most runoff flows through the ice sheet and then along the ice sheet bed to a land or ocean outlet. I have a routing code that addresses this and provides better output locations for greenland and antarctica (in the rare cases it has it) runoff. However, it is offline. @zyj8881357, can we touch base via a teams meeting to coordinate how to approach this? It might be the right time to develop catchments for the ice sheets. We have the data, I just need to do the work.

All, @zyj8881357, @biljanaorescanin, apologies for the long delay. I have produced outlets and associated drainage basins for Greenland that are gridded to the GEOS510arcsec_mask.nc. The steps closely, but not completely follow that of Verdin (2013, Final Report: High Resolution Topographic Analysis for GMAO’s Catchment LSM) and use ArcticDEM v3 100m mosaic cropped to Greenland (Porter et al., 2018; https://www.pgc.umn.edu/data/arcticdem/). 100 m is approximately 10 arc seconds at a latitude of 70° N. Steps 3-6 were completed using Topotoolbox v2 (Schwanghart & Scherler, 2014; https://topotoolbox.wordpress.com/).

1.     Crop DEM Greenland only. 2.     Fill a small number of missing grid cells using bilinear interpolation. No extrapolation was used, and this step was only designed to remove NaN values assigned to a handful of locations across Greenland. 4.     Hydrologically condition DEM by filling sinks. 5.     Calculate flow direction and flow accumulation grids using the D8 routing algorithm. 6.     Delineate ‘surface streams’ and drainage basins using an area threshold of 250 km2, the same threshold as in Verdin (2013). Streams are extracted to vector format and contain the Strahler stream order numbering system; drainage basins retain a raster format. See below for streams and 7.     Extract drainage basin outlets. For Greenland, at the given resolution there are 525 outlet locations. 8.     Map drainage basins and outlets to the GEOS10arcsec_mask.nc grid using ‘nearest’. This mapping follows the methodology outlined in Verdin (2013) & Mahanama et al. (2015; GMAO Tech Report 39). 9.     Confirm each point is appropriately assigned a drainage basin number.   Notes (in no particular order) 1.     @zyj8881357 expect an email about this and about how you would like the data transferred shortly. 2.     I have done this for all of Greenland (both ice and land) because Yujin indicated a lack of outlets across Greenland and because previous efforts excluded Greenland. I can also produce drainage basins and outlets for ice (or ice sheet) only points as well. This would reduce the number of outlet locations around Greenland by excluding many of the islands. 3.     Outlet locations are on land and still need to be moved to the nearest ocean cell. 4.     This information should be sufficient to create the routing table. While Pfaffstetter interbasins can be determined using the surface streams from step 5, including Greenland ice sheet Pfaffstetter coding system would likely cause issues (or require a rearrangement of the tile numbering) in subsequent processes because the 200000000 ice value would be replaced. My suggestion would be to append add a basin number to the Greenland ice values instead – though this is completely up for discussion because I do not know what the exact impacts would be. 5.     In the original routing scheme, is that water from land ice (200000000) locations is routed to the nearest ocean location using an expanding spiral search, so some modification of the code will be needed to exclude the Greenland ice points because my previous work modified the routing table after the spiral search. 6.     There should be future consideration of how to deal with glacier runoff routing. 7.     Because there is no dynamic ice sheet, water is theoretically routed across the surface. In reality, much of the water produced on the ice surface only flows short distances before flowing englacially to the ice sheet bed. The traditional routing scheme for subglacial water accounts for the overburden pressure of the ice and basal topography in the potential gradient calculation; however, to substantially differ from surface routing the bed slope must exceed the surface slope by ~10x, which is uncommon. 8.     Though the streams and drainage basins reach the ice sheet divide, runoff generally only occurs along the ice sheet edges.

Note: middle panel legend refers to stream order.

Fantastic! Thank you @lcandre2 can you please push your script(s) for steps: 1- 9 to github or somewhere else?

[gmao-rreichle]

Many thanks, @lcandre2!! Re. the scripts, they should definitely be added to github but we need to coordinate where exactly they should go.

[lcandre2]

Agreed @gmao-rreichle and @sanAkel. I have a bit of commenting to clean up, but I am open to where these should be located. I was thinking somewhere like @GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils But, I am open to all other suggestions too.

[sanAkel]

Agreed @gmao-rreichle and @sanAkel. I have a bit of commenting to clean up, but I am open to where these should be located. I was thinking somewhere like

@GEOSgcm_GridComp/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils

But, I am open to all other suggestions too.

Great! Thanks @lcandre2.

While I have no particular preference in following, just listing our choices (I can think of):

1. Your ⬆️ suggested path
2. This repo: https://github.com/GEOS-ESM/GEOS_Util
3. https://github.com/GEOS-ESM/GMAO_Shared/tree/main/GEOS_Shared

If above choices: 2 and 3 are to be considered, then it would a be under a new folder.

I don't think any existing code depends on this, because it's new!

For that reason, and to make @sdrabenh's life simpler, I lean towards 2 ⬆️.

[gmao-rreichle]

@lcandre2, @sanAkel: If I'm not mistaken, the code in question is a tool to create input for make_bcs. We do have similar code for generating land surface inputs to make_bcs here: GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/ I suggest placing @lcandre2's code into a new sub-directory ./route/landice/ FYI, there is a lot more such code that needs to be captured. We're trying to organize our work on make_bcs here: https://github.com/GEOS-ESM/GEOSgcm_GridComp/projects/4 The "preproc" element is captured under item 6 in the current To-Do list: "Identify ancillary input files to make_bcs that require preprocessing..." (The To-Do list is in need of updating and should change frequently.)

cc: @biljanaorescanin @weiyuan-jiang

[lcandre2]

@lcandre2, @sanAkel: If I'm not mistaken, the code in question is a tool to create input for make_bcs. We do have similar code for generating land surface inputs to make_bcs here: GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/preproc/ I suggest placing @lcandre2's code into a new sub-directory ./route/landice/ FYI, there is a lot more such code that needs to be captured. We're trying to organize our work on make_bcs here: https://github.com/GEOS-ESM/GEOSgcm_GridComp/projects/4 The "preproc" element is captured under item 6 in the current To-Do list: "Identify ancillary input files to make_bcs that require preprocessing..." (The To-Do list is in need of updating and should change frequently.)

cc: @biljanaorescanin @weiyuan-jiang

@sanAkel , @gmao-rreichle , I think Rolf's approach is the most logical to me since these scripts are directly related to the land surface. I'll clean them up and make a pull request with the additions.

[zyj8881357]

Hi all,

For this issue, I have a new version of the outlets locations. I have emailed to @yvikhlya about this, and also make an update here.

For the new version of the outlets locations files, the updates are:

1. Outlets from Greenland created by Lauren @lcandre2 are included.
2. The locations are represented by the lat lon indexes for the grid CF0180x6C_TM1440xTM1080-Pfafstetter.rst (43200x21600), as identified in this issue.
3. All the locations were moved to the ocean cells prescribed by the rst file CF0180x6C_TM1440xTM1080-Pfafstetter.rst .
4. Endorheic basins (with inland sinks) are not included in the files to ensure that all the outlet locations are at ocean cells. This issue should be treated in the future for the water balance issue.

Anyone can download them at: https://drive.google.com/drive/folders/1a2eS6zaAFRVVSpms7WqJbFrmrxEITVCD?usp=sharing

Three files are available:

1. latraster.nc: The file shows the lat indexes of the outlets for each cell in CF0180x6C_TM1440xTM1080-Pfafstetter.rst
2. lonraster.nc : The file shows the lon indexes of the outlets for each cell in CF0180x6C_TM1440xTM1080-Pfafstetter.rst
3. Outlet_latlon.43200x21600 : A outlet location files (unformatted) that could be directly read by make_bcs (specifically, https://github.com/GEOS-ESM/GEOSgcm_GridComp/blob/develop/GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mk_runofftbl.F90). This file is generated from latraster.nc and lonraster.nc through my read_riveroutlet.F90 (also uploaded in the package). I generated the file for your convenience. But if any change needs to be made for the locations, they should be made on the latraster.nc and lonraster.nc, and then get the Outlet_latlon.43200x21600 regenerated .

I hope this can help to fix the issue. In the future, improvements could be made on the file mk_runofftbl.F90 to let it be able to read the locations in degree (rather than in indexes as the current version) and adjust those locations based on different kinds of ocean grids.

[yvikhlya]

@zyj8881357

All the locations were moved to the ocean cells prescribed by the rst file CF0180x6C_TM1440xTM1080-Pfafstetter.rst

Does this mean MOM ocean cells or GEOS ocean cells? Do all *.TRN files for coupled grids need to be regenerated now?

[zyj8881357]

@zyj8881357

All the locations were moved to the ocean cells prescribed by the rst file CF0180x6C_TM1440xTM1080-Pfafstetter.rst

Does this mean MOM ocean cells or GEOS ocean cells? Do all *.TRN files for coupled grids need to be regenerated now?

I think they are in the MOM ocean cells now, but you may need to check this to see if further adjustments are needed. Yes, you need to use the file Outlet_latlon.43200x21600 to regenerate the *.TRN files by the make_bcs.

[yvikhlya]

No, this does not take into account ocean land-sea mask. Land-sea mask adjustment is better to do in mk_runofftbl.F90, which is still missing this part.

[yvikhlya]

@zyj8881357 There is a little inconsistency in both latraster.nc and lonraster.nc. In both files data field has attribute _FillValue = -9999, but the actual sentinel value is -999. This confuses netcdf tools.

[zyj8881357]

@zyj8881357 There is a little inconsistency in both latraster.nc and lonraster.nc. In both files data field has attribute _FillValue = -9999, but the actual sentinel value is -999. This confuses netcdf tools.

Sorry for the inconvenience. Yes, the actual missing value is set to -999 to be consistent with the value used in the mk_runofftbl.F90. But I forgot to change it in the netcdf attribution. Do you need me to fix that for the nc files, or you can do it by yourself?

[yvikhlya]

@zyj8881357 There is a little inconsistency in both latraster.nc and lonraster.nc. In both files data field has attribute _FillValue = -9999, but the actual sentinel value is -999. This confuses netcdf tools.

Sorry for the inconvenience. Yes, the actual missing value is set to -999 to be consistent with the value used in the mk_runofftbl.F90. But I forgot to change it in the netcdf attribution. Do you need me to fix that for the nc files, or you can do it by yourself?

I can fix these files if I need to use them. Just pointing out inconsistency in case somebody else wants to use them.

[sanAkel]

@zyj8881357 thank you! Can you upload an ascii file with: River (name) lat (deg N) lon (deg E) ABCD y_1 x_1 You should be able to upload it here, if not, please share via google drive; assuming no access to discover cluster at NCCS. Thanks again.

I uploaded it to the google drive. The link is: https://docs.google.com/spreadsheets/d/1FbQsqoYL-knSmF_qzKlOcC8Pp3kxW4Ok/edit?usp=share_link&ouid=110154394974887000797&rtpof=true&sd=true

It's a spreadsheet. I also added the mean discharge for each river. Please check whether these locations are consistent with the ocean mask in the ocean model.

Hi @zyj8881357, Would you be able to upload this and share a link? ⬆️ link does not work. If you prefer, you can email me. Thanks a lot!

[zyj8881357]

@zyj8881357 thank you! Can you upload an ascii file with: River (name) lat (deg N) lon (deg E) ABCD y_1 x_1 You should be able to upload it here, if not, please share via google drive; assuming no access to discover cluster at NCCS. Thanks again.

I uploaded it to the google drive. The link is: https://docs.google.com/spreadsheets/d/1FbQsqoYL-knSmF_qzKlOcC8Pp3kxW4Ok/edit?usp=share_link&ouid=110154394974887000797&rtpof=true&sd=true It's a spreadsheet. I also added the mean discharge for each river. Please check whether these locations are consistent with the ocean mask in the ocean model.

Hi @zyj8881357, Would you be able to upload this and share a link? ⬆️ link does not work. If you prefer, you can email me. Thanks a lot!

The issue should be solved in GEOSagcm_GridComp/GEOSphysics_GridComp/GEOSsurface_GridComp/Utils/Raster/makebcs/mk_runofftbl.F90 by introducing a subroutine outlets_to_ocean(Gridname,lons,lats,nx,ny) which moves the outlets on land to the nearest ocean point. Are you still working on this and wanting to see those locations? If so, I can run the code and make the spreadsheet for you again.

[sanAkel]

@zyj8881357 in Jul, I transferred to NOAA, EMC. I'm interested in looking at where your catchment points are, if you happen to have that spreadsheet handy and willing to share, it would help. Thanks!