SI3-10_EAP_rheology
Last edition: Wikinfo(changed_ts)? by Wikinfo(changed_by)?
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.
Summary
| Action | SI3-10_EAP_rheology |
|---|---|
| PI(S) | Stefanie Rynders |
| Digest | Add new module icedyn_rhg_eap.F90 with similar structure to the existing icedyn_rhg_evp.F90 module to include EAP rheology as an option in SI3. |
| Dependencies | |
| Branch | source:/NEMO/branches/{YEAR}/dev_r{REV}_{ACTION_NAME} |
| Previewer(s) | Andrew Coward |
| Reviewer(s) | Clement Rousset? |
| Ticket | #2325 |
Description
Current EVP rheology assumes ice is isotropic, but observations show linear kinematic features (LKF) which are sign of anisotropic behaviour.
Elastic anisotropic plastic (EAP) rheology has been developed by UCL to address this shortcoming (refs 1-3 theoretical development, refs 4-5 for implementation in CICE). It is routinely used in CICE simulations and shows more realistic LKF patterns.
The IMMERSE project has a milestone and deliverable to implement EAP rheology in SI3.
Implementation
Development cautions
- EAP used in CICE with explicit solver, hasn't been tested with implicit solver (check CICE manual)
- EAP hasn't been used with adaptive solver (Kimmritz) ' stability unknown ' not used for now
- EAP hasn't been used with Lemieux landfast ice changes ' stability unknown ' not used for now
- need 8 functions, though coding rules say it is to be avoided
- all recycled CICE code will need to be rewritten to match NEMO coding conventions.
Module list
modules to be changed:
trunk/src/ICE/icedyn_adv_pra.F90 trunk/src/ICE/icedyn_adv_umx.F90 trunk/src/ICE/icedyn_rhg.F90
new module
trunk/src/ICE/icedyn_rhg_eap.F90
Code modifications
EAP code based on Heorton et al. 2018 implementation in CICE
icedyn_rhg.F90
This is the main rheology module, contains two subroutines: ice_dyn_rhg and ice_dyn_rhg_init. Some placeholders for EAP rheology already defined
The subroutine ice_dyn_rhg does timing and conservation checks, calls the rheologies (select case pattern) and calls the writing of restart files
The subroutine ice_dyn_rhg_init reads the rheology parameters from the namelist and writes to the log
no namelist parameters needed for eap specifically, need to use nn_evp and rn_relast for the subcycling
actions
- select case parameter np_rhg_EAP to choose EAP rheology, defined on line 35 to be uncommented
- namelist parameter to be added for EAP after line 38
LOGICAL :: ln_rhg_EAP ! EAP rheology
- np_rhg_EAP used to switch on eap rheology on line 139 to be uncommented
- namelist parameter to be added in line 112
NAMELIST/namdyn_rhg/ ln_rhg_EVP, ln_aEVP, rn_creepl, rn_ecc , nn_nevp, rn_relast, ln_rhg_EAP
- line of log output to be added for eap, proposed after line 133 ?
WRITE(numout,*) ' rheology EAP (icedyn_rhg_eap) ln_rhg_EAP = ', ln_rhg_EAP
- to be added USE statement, after line 18
USE icedyn_rhg_eap ! sea-ice: EAP rheology
- to be added CASE option, after line 82
CASE ( np_rhgEAP)
CALL ice_dyn_rhg_eap( kt, stress1_i, stress2_i, stress12_i, shear_i, divu_i, delta_i )
- to be added call to restart file writing in ice_dyn_rhg, after line 86
IF( ln_rhg_EAP ) CALL rhg_eap_rst( 'WRITE', kt )
- to be added call to restart file reading in ice_dyn_rhg_init, after line 142
IF( ln_rhg_EAP ) CALL rhg_eap_rst( 'READ' ) !* read or initialize all required files
icedyn_rhg_eap.F90
new module to be added
The structure will be the same as icedyn_rhg_evp.F90, with some code and extra functions added The subroutine ice_dyn_rhg_eap contains the stress calculation and momentum balance and subcycling
- compute mask at F point: no changes needed
- compute ice, snow mass + ice strength: initialisation of variables to be added: uses 8 functions for the variables w1, w2, s11kr, s12kr, s22r, s11ks, s12ks, s22ks
- keep landfast ice option, but test later unsure whether extra tensile strength will be stable hasn't been used with EAP in other models
- compute terms of momentum equation: wind/ocean stress, mass terms, coriolis term: no changes needed
- solve momentum equation iteratively: use Bouillon option (more stable), possibly test Kimmritz option later
- insert stress update: insert structure tension evolution (stepa subroutine), line 688 unsure whether this will be stable with implicit solver instead of explicit one ->ask Harry (note in stepa says it is implicit, clarify)
- calls calc_ffrac 4 times, 4 corners of cell and average: ( warning: s11, s22 not defined on same point as s12 in SI3 need 4 corners, check?)
- calculate delta, shear, divergence for redistribution: no changes needed
- diagnostics: add output of anisotropy with iom_put calls need to define fields in xml files
- The subroutine rhg_eap_rst: add structure tensor components (anistropy) restart writing/reading variables a11_1, a11_2, a11_3, a11_4 and a12_1, a12_2, a12_3, a12_4; pay attention to order of reads 8 functions for yield surface lookup table inside functions everything scalar
- calc_ffrac subroutine? uses stresses and anisotropy tensor (both defined on corners), but inside function everything scalar - separate subroutine will not be needed if only one grid point is enough
icedyn_adv_pra.F90
The first advection module, containing tracers: add anisotropy tracer
icedyn_adv_umx.F90
The second advection module, containing tracers: add anisotropy tracer
References
1) Wilchinsky, A.V. and Feltham, D.L. A continuum anisotropic model of sea-ice dynamics. Proc. R. Soc. Lond. A, 460(2047): 2105-2140, July 2004. doi: 10.1098/rspa.2004.1282 .
2) Wilchinsky, A.V. and Feltham, D.L. Dependence of Sea Ice Yield-Curve Shape on Ice Thickness. J. Phys. Oceanogr., 34(12): 2852-2856, 2004. doi: 10.1175/JPO2667.1
3) Wilchinsky, A.V. and Feltham, D.L. Modelling the rheology of sea ice as a collection of diamond-shaped floes. J. Nonnewton. Fluid Mech., 138(1): 22-32, September 2006. doi: 10.1016/j.jnnfm.2006.05.001.
4) Tsamados, M.; Feltham, D. L., and Wilchinsky, A. V. Impact of a new anisotropic rheology on simulations of Arctic sea ice. J. Geophys. Res., 118 (1):91-107, January 2013. doi: 10.1029/2012JC007990.
5) Heorton, H. D. B. S. and Feltham, D. L. and Tsamados, M. Stress and deformation characteristics of sea ice in a high-resolution, anisotropic sea ice model. Phil. Trans. R. Soc. A, 376(2129), 2018. doi: 10.1098/rsta.2017.0349
Test case
We will use the ice_adv2d testcase modified to examine the effect of the rheology by imposing different controlled boundary conditions. The approach is inspired by the idealised simulations in Heorton et al. 2018.
Documentation updates
Section on EAP rheology to be included SI3 reference manual section 3 : ice dynamics.
Preview
update: EAP rheology has an impact on the closing rate for ridging/rafting. Currently this is hardcoded appropriate for EVP in icedyn_rdgrft.F90. An option will be added in icedyn_rdgrft.F90 to use EAP rate instead. update 2: Current implementation in CICE ignores the impact of advection on the anisotropy. Initial implementation in SI3 will do the same, so no extra tracer is needed for now. This can be reasessed later.
Tests
The test case is called ICE_RHEO, included in the development branch in the tests folder. It is based on the ice_adv2d test case, but with some differences in the forcing and boundary conditions, see the usrdef files in the MY_SRC directory. There is no single verification value, rather the intention is to examine the deformation patterns caused by different rheologies, which can be compared to those in Heorton et al. 2018 (CICE-based). The rheologies are EVP vs EAP, but it could be used for other future rheology options too.
Current status: there is still an issue preventing to reproduce Harry's results, probably in the boundary conditions, which I need to recheck with him.
Convergence of the velocities can be checked with code from changeset 12741 which is not yet incorporated in this branch.
Review
...
