nonlinear-freesurf.rst

Non-linear free-surface

Options have been added to the model that concern the free surface formulation.

Pressure/geo-potential and free surface

For the atmosphere, since ϕ = ϕ_topo − ∫_ps^pαdp, subtracting the reference state defined in section hpe-p-geo-potential-split :

$$\phi_o = \phi_{topo} - \int^p_{p_o} \alpha_o dp \hspace{5mm}\mathrm{with}\hspace{3mm} \phi_o(p_o)=\phi_{topo}$$

we get:

ϕ′ = ϕ − ϕ_o = ∫_p^p_sαdp − ∫_p^p_oα_odp

For the ocean, the reference state is simpler since ρ_c does not dependent on z (b_o = g) and the surface reference position is uniformly z = 0 (R_o = 0), and the same subtraction leads to a similar relation. For both fluids, using the isomorphic notations, we can write:

ϕ′ = ∫_r^r_surfb dr − ∫_r^R_ob_odr

and re-write as:

ϕ′ = ∫_Ro^r_surfb dr + ∫_r^R_o(b − b_o)dr

or:

ϕ′ = ∫_Ro^r_surfb_odr + ∫_r^r_surf(b − b_o)dr

In section finding_the_pressure_field, following eq. split-phi-Ro, the pressure/geo-potential ϕ′ has been separated into surface (ϕ_s), and hydrostatic anomaly (ϕ′_hyd). In this section, the split between ϕ_s and ϕ′_hyd is made according to equation split-phi-bo. This slightly different definition reflects the actual implementation in the code and is valid for both linear and non-linear free-surface formulation, in both r-coordinate and r*-coordinate.

Because the linear free-surface approximation ignores the tracer content of the fluid parcel between R_o and r_surf = R_o + η, for consistency reasons, this part is also neglected in ϕ′_hyd :

ϕ′_hyd = ∫_r^r_surf(b − b_o)dr ≃ ∫_r^R_o(b − b_o)dr

Note that in this case, the two definitions of ϕ_s and ϕ′_hyd from equations split-phi-Ro and split-phi-bo converge toward the same (approximated) expressions: ϕ_s = ∫_Ro^r_surfb_odr and ϕ′_hyd = ∫_r^R_ob′dr. On the contrary, the unapproximated formulation (see free_surf_effect_col_thick) retains the full expression: ϕ′_hyd = ∫_r^r_surf(b − b_o)dr . This is obtained by selecting nonlinFreeSurf =4 in parameter file data. Regarding the surface potential:

$$\phi_s = \int_{R_o}^{R_o+\eta} b_o dr = b_s \eta \hspace{5mm}\mathrm{with}\hspace{5mm} b_s = \frac{1}{\eta} \int_{R_o}^{R_o+\eta} b_o dr$$

b_s ≃ b_o(R_o) is an excellent approximation (better than the usual numerical truncation, since generally |η| is smaller than the vertical grid increment).

For the ocean, ϕ_s = gη and b_s = g is uniform. For the atmosphere, however, because of topographic effects, the reference surface pressure R_o = p_o has large spatial variations that are responsible for significant b_s variations (from 0.8 to 1.2 [m³/kg]). For this reason, when uniformLin_PhiSurf =.FALSE. (parameter file data, namelist PARAM01) a non-uniform linear coefficient b_s is used and computed (INI_LINEAR_PHISURF <model/src/ini_linear_phisurf.F>) according to the reference surface pressure p_o: b_s = b_o(R_o) = c_pκ(p_o/P_SL^o)^{(κ − 1)}θ_ref(p_o), with P_SL^o the mean sea-level pressure.

Free surface effect on column total thickness (Non-linear free-surface)

The total thickness of the fluid column is r_surf − R_fixed = η + R_o − R_fixed. In most applications, the free surface displacements are small compared to the total thickness η ≪ H_o = R_o − R_fixed. In the previous sections and in older version of the model, the linearized free-surface approximation was made, assuming r_surf − R_fixed ≃ H_o when computing horizontal transports, either in the continuity equation or in tracer and momentum advection terms. This approximation is dropped when using the non-linear free-surface formulation and the total thickness, including the time varying part η, is considered when computing horizontal transports. Implications for the barotropic part are presented hereafter. In section tracer-cons-nonlinear-freesurface consequences for tracer conservation is briefly discussed (more details can be found in Campin et al. (2004) cam:04) ; the general time-stepping is presented in section nonlin-freesurf-timestepping with some limitations regarding the vertical resolution in section nonlin-freesurf-dzsurf.

In the non-linear formulation, the continuous form of the model equations remains unchanged, except for the 2D continuity equation discrete-time-backward-free-surface which is now integrated from R_fixed(x, y) up to r_surf = R_o + η :

$$\epsilon_{fs} \partial_t \eta = \left. \dot{r} \right|_{r=r_{surf}} + \epsilon_{fw} ({\mathcal{P-E}}) = - {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o+\eta} \vec{\bf v} dr + \epsilon_{fw} ({\mathcal{P-E}})$$

Since η has a direct effect on the horizontal velocity (through ∇_hΦ_surf), this adds a non-linear term to the free surface equation. Several options for the time discretization of this non-linear part can be considered, as detailed below.

If the column thickness is evaluated at time step n, and with implicit treatment of the surface potential gradient, equations eq-solve2D and eq-solve2D_rhs become:

$$\begin{aligned} \epsilon_{fs} {\eta}^{n+1} - {\bf \nabla}_h \cdot \Delta t^2 (\eta^{n}+R_o-R_{fixed}) {\bf \nabla}_h b_s {\eta}^{n+1} = {\eta}^*\end{aligned}$$

where

$$\begin{aligned} {\eta}^* = \epsilon_{fs} \: {\eta}^{n} - \Delta t {\bf \nabla}_h \cdot \int_{R_{fixed}}^{R_o+\eta^n} \vec{\bf v}^* dr \: + \: \epsilon_{fw} \Delta_t ({\mathcal{P-E}})^{n}\end{aligned}$$

This method requires us to update the solver matrix at each time step.

Alternatively, the non-linear contribution can be evaluated fully explicitly:

$$\begin{aligned} \epsilon_{fs} {\eta}^{n+1} - {\bf \nabla}_h \cdot \Delta t^2 (R_o-R_{fixed}) {\bf \nabla}_h b_s {\eta}^{n+1} = {\eta}^* +{\bf \nabla}_h \cdot \Delta t^2 (\eta^{n}) {\bf \nabla}_h b_s {\eta}^{n}\end{aligned}$$

This formulation allows one to keep the initial solver matrix unchanged though throughout the integration, since the non-linear free surface only affects the RHS.

Finally, another option is a “linearized” formulation where the total column thickness appears only in the integral term of the RHS eq-solve2D_rhs but not directly in the equation eq-solve2D.

Those different options (see nonlinFreeSurf-flags) have been tested and show little differences. However, we recommend the use of the most precise method (nonlinFreeSurf =4) since the computation cost involved in the solver matrix update is negligible.

Non-linear free-surface flags

Parameter	Value	Description
nonlinFreeSurf	-1	linear free-surface, restart from a pickup file produced with #undef EXACT_CONSERV code
	0	linear free-surface (= default)
	4	full non-linear free-surface
	3	same as 4 but neglecting ∫RoRo + ηb′dr in Φ′hyd
	2	same as 3 but do not update cg2d solver matrix
	1	same as 2 but treat momentum as in linear free-surface
select_rStar	0	do not use r* vertical coordinate (= default)
	2	use r* vertical coordinate
	1	same as 2 but without the contribution of the slope of the coordinate in ∇Φ

Tracer conservation with non-linear free-surface

To ensure global tracer conservation (i.e., the total amount) as well as local conservation, the change in the surface level thickness must be consistent with the way the continuity equation is integrated, both in the barotropic part (to find η) and baroclinic part (to find w = ṙ).

To illustrate this, consider the shallow water model, with a source of fresh water (𝒫):

$$\partial_t h + \nabla \cdot h \vec{\bf v} = \mathcal{P}$$

where h is the total thickness of the water column. To conserve the tracer θ we have to discretize:

$$\partial_t (h \theta) + \nabla \cdot ( h \theta \vec{\bf v}) = \mathcal{P} \theta_{\mathrm{rain}}$$

Using the implicit (non-linear) free surface described above (press_meth_linear) we have:

$$\begin{aligned} \begin{aligned} h^{n+1} = h^{n} - \Delta t \nabla \cdot (h^n \, \vec{\bf v}^{n+1} ) + \Delta t \mathcal{P} \\\end{aligned} \end{aligned}$$

The discretized form of the tracer equation must adopt the same “form” in the computation of tracer fluxes, that is, the same value of h, as used in the continuity equation:

$$\begin{aligned} h^{n+1} \, \theta^{n+1} = h^n \, \theta^n - \Delta t \nabla \cdot (h^n \, \theta^n \, \vec{\bf v}^{n+1}) + \Delta t \mathcal{P} \theta_{rain}\end{aligned}$$

The use of a 3 time-levels time-stepping scheme such as the Adams-Bashforth make the conservation sightly tricky. The current implementation with the Adams-Bashforth time-stepping provides an exact local conservation and prevents any drift in the global tracer content (Campin et al. (2004) cam:04). Compared to the linear free-surface method, an additional step is required: the variation of the water column thickness (from hⁿ to h^n + 1) is not incorporated directly into the tracer equation. Instead, the model uses the G_θ terms (first step) as in the linear free surface formulation (with the “surface correction” turned “on”, see tracer section):

$$G_\theta^n = \left(- \nabla \cdot (h^n \, \theta^n \, \vec{\bf v}^{n+1}) - \dot{r}_{surf}^{n+1} \theta^n \right) / h^n$$

Then, in a second step, the thickness variation (expansion/reduction) is taken into account:

$$\theta^{n+1} = \theta^n + \Delta t \frac{h^n}{h^{n+1}} \left( G_\theta^{(n+1/2)} + \mathcal{P} (\theta_{\mathrm{rain}} - \theta^n )/h^n \right)$$

Note that with a simple forward time step (no Adams-Bashforth), these two formulations are equivalent, since $(h^{n+1} - h^{n})/ \Delta t = \mathcal{P} - \nabla \cdot (h^n \, \vec{\bf v}^{n+1} ) = P + \dot{r}_{surf}^{n+1}$

Time stepping implementation of the non-linear free-surface

The grid cell thickness was hold constant with the linear free-surface; with the non-linear free-surface, it is now varying in time, at least at the surface level. This implies some modifications of the general algorithm described earlier in sections adams-bashforth-sync and adams-bashforth-staggered.

A simplified version of the staggered in time, non-linear free-surface algorithm is detailed hereafter, and can be compared to the equivalent linear free-surface case (eq. Gv-n-staggered to t-n+1-staggered) and can also be easily transposed to the synchronous time-stepping case. Among the simplifications, salinity equation, implicit operator and detailed elliptic equation are omitted. Surface forcing is explicitly written as fluxes of temperature, fresh water and momentum, $\mathcal{Q}^{n+1/2}, \mathcal{P}^{n+1/2}, \mathcal{F}_{\bf v}^n$ respectively. hⁿ and dhⁿ are the column and grid box thickness in r-coordinate.

ϕ_hydⁿ = ∫b(θⁿ, Sⁿ, r)dr

$$\vec{\bf G}_{\vec{\bf v}}^{n-1/2}\hspace{-2mm} = \vec{\bf G}_{\vec{\bf v}} (dh^{n-1},\vec{\bf v}^{n-1/2}) \hspace{+2mm};\hspace{+2mm} \vec{\bf G}_{\vec{\bf v}}^{(n)} = \frac{3}{2} \vec{\bf G}_{\vec{\bf v}}^{n-1/2} - \frac{1}{2} \vec{\bf G}_{\vec{\bf v}}^{n-3/2}$$

$$\vec{\bf v}^{*} = \vec{\bf v}^{n-1/2} + \Delta t \frac{dh^{n-1}}{dh^{n}} \left( \vec{\bf G}_{\vec{\bf v}}^{(n)} + F_{\vec{\bf v}}^{n}/dh^{n-1} \right) - \Delta t \nabla \phi_{hyd}^{n}$$

$$\longrightarrow update \phantom{x} model \phantom{x} geometry : {\bf hFac}(dh^n)$$

$$\begin{aligned} \begin{aligned} \eta^{n+1/2} \hspace{-1mm} & = \eta^{n-1/2} + \Delta t P^{n+1/2} - \Delta t \nabla \cdot \int \vec{\bf v}^{n+1/2} dh^{n} \\\ & = \eta^{n-1/2} + \Delta t P^{n+1/2} - \Delta t \nabla \cdot \int \!\!\! \left( \vec{\bf v}^* - g \Delta t \nabla \eta^{n+1/2} \right) dh^{n}\end{aligned} \end{aligned}$$

$$\vec{\bf v}^{n+1/2}\hspace{-2mm} = \vec{\bf v}^{*} - g \Delta t \nabla \eta^{n+1/2}$$

$$h^{n+1} = h^{n} + \Delta t P^{n+1/2} - \Delta t \nabla \cdot \int \vec{\bf v}^{n+1/2} dh^{n}$$

$$G_{\theta}^{n} = G_{\theta} ( dh^{n}, u^{n+1/2}, \theta^{n} ) \hspace{+2mm};\hspace{+2mm} G_{\theta}^{(n+1/2)} = \frac{3}{2} G_{\theta}^{n} - \frac{1}{2} G_{\theta}^{n-1}$$

$$\theta^{n+1} =\theta^{n} + \Delta t \frac{dh^n}{dh^{n+1}} \left( G_{\theta}^{(n+1/2)} +( P^{n+1/2} (\theta_{\mathrm{rain}}-\theta^n) + \mathcal{Q}^{n+1/2})/dh^n \right) \nonumber$$

Two steps have been added to linear free-surface algorithm (eq. Gv-n-staggered to t-n+1-staggered): Firstly, the model “geometry” (here the hFacC,W,S) is updated just before entering SOLVE_FOR_PRESSURE <model/src/solve_for_pressure.F>, using the current dhⁿ field. Secondly, the vertically integrated continuity equation h-n+1-nlfs has been added (exactConserv =.TRUE., in parameter file data, namelist PARM01) just before computing the vertical velocity, in subroutine INTEGR_CONTINUITY <model/src/integr_continuity.F>. Although this equation might appear redundant with nstar-nlfs, the integrated column thickness h^n + 1 will be different from η^n + 1/2 + H  in the following cases:

• when Crank-Nicolson time-stepping is used (see crank-nicolson_baro).
• when filters are applied to the flow field, after v-n+1-nlfs, and alter the divergence of the flow.
• when the solver does not iterate until convergence; for example, because a too large residual target was set (cg2dTargetResidual, parameter file data, namelist PARM02).

In this staggered time-stepping algorithm, the momentum tendencies are computed using dh^n − 1 geometry factors Gv-n-nlfs and then rescaled in subroutine TIMESTEP <model/src/timestep.F>, vstar-nlfs, similarly to tracer tendencies (see tracer-cons-nonlinear-freesurface). The tracers are stepped forward later, using the recently updated flow field ${\bf v}^{n+1/2}$ and the corresponding model geometry dhⁿ to compute the tendencies Gt-n-nlfs; then the tendencies are rescaled by dhⁿ/dh^n + 1 to derive the new tracers values (θ, S)^n + 1 (t-n+1-nlfs, in subroutines CALC_GT <model/src/calc_gt.F>, CALC_GS <model/src/calc_gs.F>).

Note that the fresh-water input is added in a consistent way in the continuity equation and in the tracer equation, taking into account the fresh-water temperature θ_rain.

Regarding the restart procedure, two 2D fields h^n − 1 and (hⁿ − h^n − 1)/Δt in addition to the standard state variables and tendencies (η^{n − 1/2}, ${\bf v}^{n-1/2}$, θⁿ, Sⁿ, ${\bf G}_{\bf v}^{n-3/2}$, G_θ, S^n − 1) are stored in a “pickup” file. The model restarts reading this pickup file, then updates the model geometry according to h^n − 1, and compute hⁿ and the vertical velocity before starting the main calling sequence (eq. phi-hyd-nlfs to t-n+1-nlfs, FORWARD_STEP <model/src/forward_step.F>).

S/R INTEGR_CONTINUITY <model/src/integr_continuity.F>

h^n + 1 − H_o : etaH ( DYNVARS.h <model/inc/DYNVARS.h> )
hⁿ − H_o : etaHnm1 ( SURFACE.h <model/inc/SURFACE.h> )
(h^n + 1 − hⁿ)/Δt : dEtaHdt ( SURFACE.h <model/inc/SURFACE.h> )

Non-linear free-surface and vertical resolution

When the amplitude of the free-surface variations becomes as large as the vertical resolution near the surface, the surface layer thickness can decrease to nearly zero or can even vanish completely. This later possibility has not been implemented, and a minimum relative thickness is imposed (hFacInf, parameter file data, namelist PARM01) to prevent numerical instabilities caused by very thin surface level.

A better alternative to the vanishing level problem relies on a different vertical coordinate r^* : The time variation of the total column thickness becomes part of the r^* coordinate motion, as in a σ_z, σ_p model, but the fixed part related to topography is treated as in a height or pressure coordinate model. A complete description is given in Adcroft and Campin (2004) adcroft:04a.

The time-stepping implementation of the r^* coordinate is identical to the non-linear free-surface in r coordinate, and differences appear only in the spacial discretization.