Wetting and Drying development work

This is the color code for the fulfilment of this form:

PI(S)

Previewer(s)

Reviewer(s)


The PI is responsible to closely follow the progress of the action, and especially to contact NEMO project manager if the delay on preview (or review) are longer than the 2 weeks expected.

Abstract

A Wetting and Drying (W/D) Method for NEMO-shelf Based on Flux Limiter

Hedong Liu and Jason Holt (NOC, UK)

NEMO discretize its computational domain with the structured C-grid, hence its W/D algorithm needs to deal with two types of grid cells: velocity grid cells and tracer grid cells. The essences of current method are to keep the positivity of watermass and avoid unphysical overshooting of velocity.At the same time, this method must respect the properties of local and global mass conservation. There are quite many existing methods for W/D,(see Balzano,A, 1998 and Medeiros,S.C. and Hagen,S.C, 2013). But there are still some robustness issues, such as model stability, computational efficiency. In the current method, a flux limiter has been developed to take care of the mass positivity while a new barotropic pressure gradient filter is developed to prevent the velocity overshooting at the wetting and drying cells on a steep slope.

In the following parts of this note, we will first derive the flux limiter and demonstrate its applications to the continuity equation, tracer equations. As the momentum equations is also linked with the W/D processes, we then detail the design of the gravity forcing filter which aims to remove the errors generated by the normal pressure gradient calculations in the W/D area. The implementation of the method in NEMO v3.6_STABLE source code will be illustrated at the last.

Flux limiter

In NEMO, the water column can be considered as a control volume V. The horizontal advection is the main process causing the wetting and drying on the grid cells (with the evaporation and precipitation as the other two reasons which will be dealt separately). For such a control volume V, we can define a flux limiter, maths fig 1, as follows:

The overshooting flux across a W/D grid cell is:

maths fig 2

To preserve the positivity of water mass at each cell, the limiter coefficient maths fig 1 must satisfy:

flux limiter is needed if:

So:

is the time step, U is the vertical integrated velocity, q is the source rate from precipitation, evaporation, and river inflow, Fi is water flux rate across each boundary of the control volume, which is positive for outflow and negative for inflow. The subscript i identifies each boundary  in the west, east, north, and south directions. Dmin is the minimal water depth to be kept in a dry grid cell while D is the total depth of the grid cell. is the horizontal area of the water column. sign() is the sign function. Obviously, should have the following property:

In other words, the strength of limiter varies from 0 to 1 which corresponding to the fully shutting down the outflow and no limiter on outflow respectively.

One case needs to be paid attention here is the above equations can produce a negative value of , This is corresponding to a large evaporation rate. In that case, we need to turn off the evaporation on that grid cell during that time step.

The limited flux (only in the horizontal directions) has the following form:

The above equation is applied to all the continuity equation and tracer equations. The resultant water column on each grid cell will have its depth greater than or equals to the predefined minimal depth Dmin.

Continuity Equation and Tracer Equations.

We need to be careful with the explicit time splitting method when the flux limiter is applied to the continuity equation.There are many external time steps within each internal timestep.The multi- barotropic stepping updates the sea surface height and vertical integrated horizontal velocity. It is tedious to keep the consistency between these two modes on the W/D grid cells. The reason to do the explicit time splitting is because of the big ratio of gravity wave speed to the baroclinic velocity. We argue here that this ratio must be a small value (less than or equals to 1) in the W/D area because of the relatively small water depth. So, we simplify this case by turning off the time splitting procedure on the W/D grid cells.

A few iterations may be needed to take into account of the interactions between the two neighbouring cells when the limiter coefficient is calculated. As this flux limiter only effects around the W/D front, the iteration should converge very quickly. Also, the iteration is only needed for the continuity equation and should not be a very much extra burden to the whole computation.

For the tracer Equation, it can simply be implemented by applying the flux limiter coefficient to the horizontal velocities. In such way, we don’t need change anything in the original source code except update the velocities.

Momentum equations

For the momentum equations, the grid cell water depth is just the weight averaged depth of the two neighboring water depths, hence it always keeps positive as long as the water depth on the tracer grid cells is positive. So, we don’t need to apply any momentum flux limiter.

The critical issue is how to maintain the model stability. The velocity nodes can be categorized into ocean nodes, land nodes, and W/D nodes based on whether wetting/drying happens or not on it. For the ocean nodes, which never got dried out, we don’t need to implement any W/D work there. For the land nodes, which never got flooded with water, we can just set the velocity to zero and don’t need to solve any equations. For the W/D nodes, some specific improvements need to be implemented to keep model stable. These improvements include filtering out the unrealistic hydrostatic pressure gradient and other nonlinear and the Coriolis forcing term.

A momentum grid cell which is between a wet grid cell and a dry grid cell, will have modified barotropic pressure gradients as demonstrated in the following two equations:

where:

Where, sign() is the Fortran sign function. D is total depth, is the sea surface height. g is gravity acceleration. P is the barotropic pressure gradient.

The bottom friction on DDW will automatically increase when the “log layer” configuration is used. The vertical viscous/diffusive coefficient on DDW will be increased manually to keep the model stable.

The surface wind stress and surface heating/cooling will be neglected at DDW cells. The justification is that these cells are thought fully or nearly dried.

Code implementation:

One problem with to implement the current method is the sigma coordinate redistribution with time. The existing NEMO code uses coefficient “mut(uvf)(:,:,:)” which is based on the physical initial vertical length scale factors. In wetting/drying case, people might encounter zero depth in some grid cells at the initiation stage. So, it is modified to use the initial sigma coordinate vertical scale to do the redistribution.

The current method tries to be independent of the existing code as much as possible, that is to say, to minimise the interferences from the wetting/drying module to the existing code. CPP method seems to be the most ideal one. But we were reminded that NEMO development policy is trying to reduce the use of CPP. So, the code is implemented with some extra coefficients which have value of 1 or zero when applied to the non-wetting/drying case.In some other parts, the wetting/drying code was added in some IF or CASE blocks. We are currently working on the implementation of this method onto NEMO v3.6 r:6397. Will provide a detailed list of the code modifications when the implementation is finished.

This section should be completed before starting to develop the code, in order to find agreement on the method beforehand.

Details

Action

PI(S)

Hedong Liu

Ticket

#1701

Branch

branches/2016/dev_r6393_NOC_WAD

Previewer(s)

Rachel Furner

Reviewer(s)

Rachel Furner

To enabling the !ticket and the source links related to your action, edit the form like a ordinary wiki page to hardcode them inside the table

Description

Implementation

Reference manual updates

Updated on 08/14/2020 19:20:56 by anonymous

Once the PI has completed this section, he should send a mail to the previewer(s) asking them to preview the work within two weeks.

Preview

Since the preview step must be completed before the PI starts the coding, the previewer(s) answers are expected to be completed within the two weeks after the PI has sent his request.
For each question, an iterative process should take place between PI and previewer(s) in order to reach a "YES" answer for each of the following questions.

Questions Answer Comment
Does the previewer agree with the proposed methodology?
Does the previewer agree with the proposed flowchart and list of routines to be changed?
Does the previewer agree with the proposed new list of variables, including agreement with coding rules?
Does the previewer agree with the proposed summary of updates in reference manual?
… … …

Updated on 08/14/2020 19:20:56 by anonymous

Once all "YES" have been reached, the PI can start the development into his development branch.

Tests

Once the development is done, the PI should complete this section below and ask the reviewers to start their review in the lower section.

Questions Answer Comment
Can this change be shown to produce expected impact? (if option activated)?
Can this change be shown to have a null impact? (if option not activated)
Detailed results of restartability and reproducibility when the option is activated. Please indicate the configuration used for this test
Detailed results of SETTE tests (restartability and reproducibility for each of the reference configuration)
Results of the required bit comparability tests been run: Are there no differences when activating the development?
If some differences appear, is reason for the change valid/understood?
If some differences appear, is the !ticket describing in detail the impact this change will have on model configurations?
Is this change expected to preserve all diagnostics?
If no, is reason for the change valid/understood?
Are there significant changes in run time/memory?
… … …

Updated on 08/14/2020 19:20:56 by anonymous

Review

A successful review is needed to schedule the merge of this development into the future NEMO release during next Merge Party (usually in November).

Code changes and documentation

Question Answer Comment
Is the proposed methodology now implemented?
Are the code changes in agreement with the flowchart defined at Preview step?
Are the code changes in agreement with list of routines and variables as proposed at Preview step?
If not, are the discrepancies acceptable?
Is the in-line documentation accurate and sufficient?
Do the code changes comply with NEMO coding standards?
Is the !ticket of development documented with sufficient details for others to understand the impact of the change?
Are the reference manual tex files now updated following the proposed summary in preview section?
Is there a need for some documentation on the web pages (in addition to in-line and reference manual)?
If yes, please describe and ask PI. A yes answer must include all documentation available.
… … …

Review Summary
Is the review fully successful?

Updated on 08/14/2020 19:20:56 by anonymous

Once review is successful, the development must be scheduled for merge during next Merge Party Meeting.

Notes from meetings

8 June 2016

Attending: Mike Bell, Andrew Coward, Jason Holt, Hedong Liu

  1. Main points
  1. History and Objectives: Jason and Hedong have previous experience of WAD in POLCOMS and FVCOM respectively. They have been working on WAD in NEMO since 2012. WAD is a key capability for NOC and it would be valuable for NOC to have it in NEMO. The Met Office needs WAD to work in NEMO for an upgrade to CMEMS north-west European shelf system expected to be required within 1-2 years. The aim is to implement a non intrusive, accurate, stable and efficient formulation of WAD. In order to be non intrusive a scheme using a minimum depth of water (@ 5-10 cm) at all potentially wet points is being implemented (similar to POM and ROMS implementation). We should publish results of our implementation. To do this we will need to compare results with other schemes (e.g. POM and ROMS implementations) for simple test cases.
  1. Review of WAD code: The main changes to an ocean model and issues associated with an implementation of WAD are:
  • a time-varying land-sea mask (wdmask) is required; the initial set up of the grid needs to allow negative depths (domgr)
  • the depth of water at ocean points cannot be allowed to go negative; some care with flux limiters is needed to ensure this. The splitting of the barotropic and baroclinic timesteps is complicated by cells wetting and drying within the baroclinic time-step
  • the thin film of water at dry cells results in an artificial horizontal pressure gradient (hpg) along a sloping bottom
  • some implementations (e.g. POM) have adjusted the bottom friction or anticipated problems with (uncoupled) surface heat fluxes generating "boiling" water of 5 cm depth
  • the 3D equations can be unstable in a very thin layer of water (with many levels); e.g. problems with CFL for vertical advection or diffusion

  1. The WAD implementation by ROMS is particularly simple. The flux limiter simply specifies that there is no flux out of a dry cell. This allows

wetting of dry cells and (in the no wind case) ensures that the artificial hpgs have no impact on the model

  1. The time-step splitting could give stability problems. In order to separate issues it would be valuable to run without time-step splitting.

We can do this by removing the updates calculated from the barotropic time-steps that are made on the baroclinic time-step.

  1. Andrew presented results using Hedong's code for a tidally driven sloping channel. It seems that when the tide is going out, strong horizontal

sea surface height gradients develop (before the layers dry out). Either the pressure gradients are being mis-calculated or the outward fluxes of water are being mis-calculated.

  1. Actions
  1. Jason to seek agreement for Hedong to focus on W&D over the next few weeks
  1. Hedong to circulate links to the papers documenting the ROMS and POM WAD implementation
  1. Hedong to
  • transfer Andrew's tidally driven sloping channel test case on his system and reproduce Andrew's results
  • run with more output diagnostics (e.g. fluxes at faces) and diagnose the cause of the problem discussed at the meeting
  • compare the parameter settings Andrew used for the test case with those in the ROMS and POM papers and adjust them to match the previous work
  • implement the ROMS flux limiter as an option within his code and use it with the test case
  • turn off the barotropic update to the baroclinic updates and re-run the test case
  • run the test case with no bottom friction (as a sensitivity test)
  1. Andrew to send instructions to Mike and Jason on how to extract the WAD code from the NEMO source so they can review its formulation
  1. Mike and Jason to
  • produce a summary of the ROMS and POM WAD algorithms
  • read Hedong's code and write a first draft documentation for it
  1. Mike to write a summary of the meeting and do a Doodle poll for a phone meeting in about 3 weeks time
Last modified 4 years ago Last modified on 2016-06-09T10:42:06+02:00

Attachments (17)

Download all attachments as: .zip