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 NEMOshelf Based on Flux Limiter
Hedong Liu and Jason Holt (NOC, UK)
NEMO discretize its computational domain with the structured Cgrid, 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, , as follows:
The overshooting flux across a W/D grid cell is:
To preserve the positivity of water mass at each cell, the limiter coefficient 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 nonwetting/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.
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.
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.
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).
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
 Main points
 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 northwest European shelf system expected to be required within 12 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 (@ 510 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.
 Review of WAD code: The main changes to an ocean model and issues associated with an implementation of WAD are:
 a timevarying landsea 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 timestep
 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
 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
 The timestep splitting could give stability problems. In order to separate issues it would be valuable to run without timestep splitting.
We can do this by removing the updates calculated from the barotropic timesteps that are made on the baroclinic timestep.
 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 miscalculated or the outward fluxes of water are being miscalculated.
 Actions
 Jason to seek agreement for Hedong to focus on W&D over the next few weeks
 Hedong to circulate links to the papers documenting the ROMS and POM WAD implementation
 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 rerun the test case
 run the test case with no bottom friction (as a sensitivity test)
 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
 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
 Mike to write a summary of the meeting and do a Doodle poll for a phone meeting in about 3 weeks time
Attachments (17)

fig01.png
(281 bytes) 
added by acc 5 years ago.
maths fig 1
 fig03.png (281 bytes)  added by acc 5 years ago.
 fig04.png (1.2 KB)  added by acc 5 years ago.
 fig06.png (2.3 KB)  added by acc 5 years ago.
 fig07.png (270 bytes)  added by acc 5 years ago.
 fig08.png (302 bytes)  added by acc 5 years ago.
 fig09.png (281 bytes)  added by acc 5 years ago.
 fig10.png (459 bytes)  added by acc 5 years ago.
 fig11.png (281 bytes)  added by acc 5 years ago.
 fig12.png (1.0 KB)  added by acc 5 years ago.
 fig13.png (723 bytes)  added by acc 5 years ago.
 fig14.png (742 bytes)  added by acc 5 years ago.
 fig15.png (1.6 KB)  added by acc 5 years ago.
 fig16.png (1.6 KB)  added by acc 5 years ago.
 fig17.png (217 bytes)  added by acc 5 years ago.
 fig05.png (818 bytes)  added by acc 5 years ago.

fig02.png
(1.3 KB) 
added by acc 5 years ago.
maths fig 2
Download all attachments as: .zip