Version 8 (modified by mike.bell, 8 years ago) (diff) 

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
Attachments (17)

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

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