On the numerical stability of surface-atmosphere coupling in weather and climate models

. Coupling the atmosphere with the underlying surface presents numerical stability challenges in cost-effective model integrations used for operational weather prediction or climate simulations. These are due to the choice of large integration time-step compared to the physical time scale of the problem, aiming at reducing computational burden, and to an explicit ﬂux coupling formulation, often preferred for its simplicity and modularity. Atmospheric models therefore use the surface-layer temperatures (representative of the uppermost soil, snow, ice, water, etc.,) at the previous integration time-step in all surface-5 atmosphere heat-ﬂux calculations and prescribe ﬂuxes to be used in the surface model integrations. Although both models may use implicit formulations for the time stepping, the explicit ﬂux coupling can still lead to instabilities. In this study, idealized simulations with a fully coupled implicit system are performed to derive an empirical relation between surface heat ﬂux and surface temperature at the new time level. Such a relation mimics the fully implicit formulation by allowing to estimate the surface temperature at the new time level without solving the surface heat diffusion problem. It is 10 based on similarity reasoning and applies to any medium with constant heat diffusion and heat capacity parameters. The advantage is that modularity of the code is maintained and that the heat ﬂux can be computed in the atmospheric model in such a way that instabilities in the snow or ice code are avoided. Applicability to snow/ice/soil models with variable density is discussed, and the loss of accuracy turns out to be small. A formal stability analysis conﬁrms that the parametrized implicit ﬂux coupling is unconditionally stable.


Introduction
Coupling atmospheric models to the underlying surface model involves both scientific and technical issues.Models of the atmospheric circulation tend to be computer intensive and therefore often employ long time steps (up to one hour), which is a challenge for stability and accuracy (Beljaars et al., 2004;Lemarié et al., 2015).The turbulent diffusion part of these codes provides the coupling to the surface, has short physical time scales near the surface and therefore needs implicit numerics for stability.The surface may be vegetation, soil, snow, ice, or a combination of these in a tile scheme.Best et al. (2004) propose a coupling strategy to the surface that has a clean interface between atmosphere and surface code, and allows to include the surface or the top part of the surface in the implicit computations.This is often necessary for stability if the physical time scale of e.g.vegetation, soil, snow or ice surface is short compared to the model time step.
The ideal solution for stability is to combine the boundary layer heat diffusion and e.g. the snow or ice layer diffusion in a single implicit solver.This has been demonstrated in a series of papers describing developments in the ORCHIDEE model (Polcher et al., 1998;Ryder et al., 2016;Wang et al., 2013).The same method was used by Schulz et al. (2001) in the Hamburg model.However, not all models do have sufficient modularity of the code to make it practical.The complication of processes like phase changes and water percolation also require implementation of conserved variables to support full implicitness (Wang et al., 2013).The standard solution is to compute fluxes at the surface on the basis of the old time level surface temperature.
It is often called "explicit flux coupling".To improve stability and accuracy West et al. (2016) recently proposed to move the flux coupling level one level down i.e. just below the surface.This has the advantage of including the fast responding surface layer in the fully implicit computations, which is beneficial for stability and accuracy.
Ongoing work at ECMWF on snow modelling raised similar issues.The existing single layer snow model (see e.g.Dutra et al., 2010), has already a minor stability issue when the snow layer becomes very thin, e.g. during the first snowfall in the season and at the final melt.This was addressed by introducing some empirical implicitness in the coupling by making an educated guess of the future snow temperature.Initial experimentation with a multilayer snow model (Dutra et al., 2012) showed even more frequent instabilities, so more implicitness in the coupling is required for stability.
In this paper, we propose a solution, that has the simplicity and modularity of the explicit flux coupling, but still has the stability of the fully implicit system.To derive simple solutions, the fully implicit coupled system is used as a reference.It is shown that the tri-diagonal set of equations corresponding to the discretized diffusion equation (for snow, ice or soil) can be converted to a relation between temperature and heat flux at the surface.The coefficients in this relation are then parametrized depending on properties of the medium, time step and vertical discretization.The coefficients are put in dimensionless form, which makes the empirical coefficients universal and applicable to any medium and any discretization.
The experimental environment in this paper is a simple model of a near surface air layer coupled to a snow pack by turbulent exchange.The atmosphere (e.g. at a height of 10 m, typical for atmospheric models) is assumed to have a diurnal cycle, and the response of temperature in the snow pack is considered.Although the following sections refer to snow only, the dimensionless framework ensures that the outcome is valid for any medium.
The following two sections (2 and 3) describe the equations for the discretized snow layer and the turbulent coupling between atmosphere and snow.Sections 4, 5 and 6 describe the numerical solution for an idealized diurnal cycle, the parametrization of the coefficients that relate heat flux and top layer snow temperature and the testing of the proposed scheme.Finally, the results and their applicability are briefly discussed in the concluding section.Also the implications of non-uniform snow density are discussed.The numerical solver and a formal stability analysis are described in Appendix A and B respectively.
We consider the diffusion equation for temperature in snow where ρ (kgm is the diffusion coefficient for heat.The boundary conditions are: For numerical stability with long time steps it is necessary to use an implicit scheme.With a vertical grid defined as in Fig. 1, the equation can be discretized as follows with boundary conditions This set of equations forms a tri-diagonal system, with diagonals A, B and C (the coefficients are defined in Appendix A).The matrix equations can be solved by successive elimination from the bottom upward such that the C-coefficients are replaced by zeros.At the same time, the equations are scaled to obtain B-coefficients that are equal to 1. Arriving at the top, it provides a solution for T n+1

1
. The solution for the other layers can be found by successive back-substitution of the temperatures going from top to bottom (see Appendix A for more details).
In case G 0 is not known, the elimination provides a linear relation between G 0 and T n+1 1 This relation can be used to achieve fully implicit coupling with the air/surface interaction formulation.

Coupling to the lowest model level of the atmosphere
To focus on stability of the atmosphere surface coupling, it is assumed that the evolution of the near atmospheric temperature is known, e.g. as in standalone simulations of the land surface.However, this is not a limitation in full 3D models that typically use an implicit solver for the turbulent diffusion.In that case the atmospheric model will perform the downward elimination process (the same way as described in Appendix B).The result is a linear relation between the n + 1 temperature at the lowest atmospheric level and the surface heat flux, which can be used with the air/land interaction formulae described below to achieve fully implicit coupling.
With a prescribed air temperature, the heat flux into the snow layer can be related to the air / surface temperature difference in the following way where G 0 is the heat flux into the snow pack, ρ a is air density, c p is air heat capacity, C H is the transfer coefficient between the atmospheric level and the surface, |U | is absolute wind speed, T a is air temperature, and T sk is temperature of the snow surface (skin temperature).
The coupling through a transfer coefficient is standard and represents the integral profile function according to Monin Obukhov (MO) similarity (see e.g.Brutsaert, 1982).The transfer coefficient in neutral conditions is related to the height of the atmospheric level, and the surface roughness lengths of momentum and heat where κ is the VonKarman constant (0.4), z a is the height of the atmospheric level, z om is the surface roughness length for momentum, and z oh is the surface roughness length for heat.Stability can be included by extending the logarithmic terms with the integral MO stability functions.
In the vertically discretized snow (see Fig. 1), the temperature of layer 1 is assumed to be at the midpoint which is different from the skin temperature.Therefore, the total conductivity between the atmosphere and the first snow layer (λ t ) is composed of two components: the turbulent transfer in the air above the surface (λ a ) and the conductivity of half of the top snow layer (λ sk ).The two conductivities are in parallel, because the inverse of conductivities (resistances) are in series, leading to the following formulation for the heat flux into the snow with Two different time stepping procedures are considered: i. Explicit flux coupling.This is the traditional approach where the expression for the surface flux uses the previous time level of the surface temperature leading to the following discretization of equation ( 11) With the explicit specification of the flux at the surface flux, the tridiagonal system can be solved directly.
ii. Implicit flux coupling.The discretization of equation ( 11) reads With this fully implicit formulation, the surface heat flux can not be specified explicitly, so it has to be found as part of the coupled atmosphere/surface system.For that purpose the tri-diagonal problem is solved in two steps.First, the elimination part is performed resulting in a solution for α and β in equation ( 8).Together with equation ( 13), T n+1 1 and G 0 can be computed: Finally the entire temperature profile can be resolved by performing the back-substitution in the tri-diagonal solver.
4 Solutions with a simple multilayer snow model In this section, solutions are considered for a 1 m thick snow layer with constant heat capacity and heat diffusion coefficients.
Idealized temperature forcing from the atmosphere is prescribed as a sinusoidal diurnal cycle.The choice of constants is documented in Table 1.The initial temperature profile at t = 0 is set to −5 o C, and a single sinusoidal diurnal cycle with an amplitude of 1 o C is imposed at the 10 m level in the atmosphere ) .
The simulations are performed with different uniform vertical discretizations and different time steps.Fig. 2 shows time series of the snow skin temperature (left column) and the ground heat flux (right column), with the two schemes.The fairly long time step of 3600 seconds is selected to illustrate stability and time truncation issues, and a short time step of 100 seconds for comparison.In the latter case time truncation errors are small for both schemes (convergence was verified).The three rows in The first thing to note is that amplitude and phase of the skin temperature diurnal cycle only have a small dependence on vertical resolution.This is surprising because the amplitude of diurnal cycle of layer 1 with ∆z = 0.2 m is only 20% of the amplitude with ∆z = 0.02 m.The reason that the skin temperature is still reasonable is due to the conductivity between the middle of the layer and the top (much lower with ∆z = 0.2 m than with ∆z = 0.02 m).So at low vertical resolution, a substantial part of the temperature signal at the snow skin is due to the "interpolation" between air and middle of the first snow layer making use of the air conductivity (λ a ) and the snow conductivity of half the top layer (λ sk ).One might interpret this result as a justification for rather low vertical resolution.However, it should be realized that the forcing has the diurnal time scale only.With faster time scales e.g.due to moving clouds and frontal passages, a relatively thick near surface layer will not be able to respond.
Although it is impossible to draw general conclusions about accuracy from limited experimentation, we note that the fully implicit solution with ∆t = 3600 s is very close to the short time step solution with ∆t = 100 s, so the long time step does not compromise accuracy in this case, although the time stepping is first order accurate only.However, the solution with explicit coupling deviates visibly from the implicit and very short time step solutions (compare the red solid curve in middle/left panel of Fig. 2 with the blue curve).Apparently, it is the mismatch of time levels in the flux computation that is detrimental to accuracy.The error is particularly visible as a phase error.
Finally, the explicit coupling turns out to be unstable for very thin snow layers (see lower panels in Fig. 2 for ∆z = 0.002.
Also for this case the long time step solution with implicit coupling is fairly accurate as it is very close to the short time step solution.These experimental results are confirmed by a formal stability analysis in Appendix B. The explicit flux coupling is unstable for a particular parameter range and the implicit flux coupling is unconditionally stable.
Because of the good stability and accuracy characteristics, we develop in the next section a parametric form of α and β in equation ( 8).
5 Scaling relations for α and β As suggested above, it is desirable to have all the flux formulations (also for the atmosphere/surface exchange) at the new time level n + 1.This implies the fully implicit option as suggested by (Polcher et al., 1998) and described in sections 2 and 3.It also requires to perform the elimination part of the tri-diagonal solver to find the relation between T n+1 1 and G 0 according to equation ( 8).For code technical reasons it is often desirable to compute the heat flux into the snow, before the snow code is actually executed.Therefore, an educated guess is made of the coefficients α and β in equation ( 8) without solving the tri-diagonal system, i.e. α and β are parametrized.
For that purpose, we make use of similarity theory for the diffusion equation with constant coefficients.If we think of an infinite medium (thick snow layer) with uniform temperature T o and make a jump at the surface to T new at t = 0, we have to consider the following basic variables: the temperature change T −T 0 at time t, T new −T 0 , K/(ρC), and depth z.According to the Buckingham Pi Theorem (Stull, 1988), 5 variables with 3 dimensions (m, s, and K), lead to two independent dimensionless groups: (T − T 0 )/(T new − T 0 ) and z/δ, where Length scale δ is the natural length scale of the medium for time scale t after which the temperature change at the surface was applied.From the physical point of view, δ is the typical depth to which the perturbation of the surface temperature has propagated at time t.The implication is that (T − T 0 )/(T new − T 0 ) is a universal function of z/δ.At this stage we do not care about the form, although the solution can be found easily by transforming the equation to the new coordinate z/δ, which allows to separate the time dependence and the depth dependence leading to an ordinary differential equations which can be solved analytically (Carslaw and Jaeger, 1959).
Similarly, we can apply an external forcing by suddenly applying a heat flux G 0 at time 0 and look for the temperature response.Instead of scaling the temperature with the temperature jump, we make the temperature change dimensionless with G 0 and obtain where h is a universal function.For z = 0, equation ( 18) is of the form of equation ( 8).With time scale ∆t and substitution of the expression for δ, we therefore expect the following scaling behavior for α It indicates the surface temperature response to a 1 W/m 2 heat flux forcing over a finite time step ∆t.
The scaling arguments above apply to the continuous system.For the discretized system, the scaling behavior of α also depends on ∆z, which introduces a dependence on the dimensionless variable ∆z/δ.For a very fine grid (∆z << δ), the discrete system behaves like the continuous system and equation ( 19) applies.For a very thick top layer (∆z >> δ), the heat flux is simply distributed over the top layer and the following applies In general the dimensionless α should be a universal function of δ/∆z, i.e.
The empirical function can be "measured" by running the numerical model as in the previous section for a range of time steps and vertical discretizations.Note that α remains constant during the time stepping and does not depend on the temperature profile.It is just a property of the tri-diagonal matrix which only contains properties of the medium, the time step and the level thickness.The results are shown in Fig. 3. Time steps range from 100 s to 3600 s, and layer thicknesses are used from 0.002 m to 0.2 m, with a total snow depth of 1 m for all simulations For small ratios of δ/∆z, the universal function should scale with equation ( 20) and for large values with (19).Surprisingly, coefficient h 0 turns out to be 1.An empirical fit is proposed that makes a smooth transition between the two regimes according to (see Fig. 3) The exponent of 1.3 has been optimized to obtain a reasonable representation of the numerical data in the transition regime.
The second parameter for which an empirical formulation is needed is β.The physical meaning of β is clear from equation (8): it is the temperature of the top snow layer at the new time level T n+1 1 in case of zero heat flux.A simple approximation would be to select the temperature of the previous time level, but this is only valid for a uniform temperature profile.For a non-uniform temperature profile, heat diffusion will homogenize temperature, which will make β different from T n 1 at the old time level.Following the scaling arguments above, we know that information propagates vertically over a distance δ during time step ∆.Therefore, we conjecture that the temperature of the old profile at depth δ is a better approximation for β than the temperature at level 1, i.e.T n δ is better than T n 1 .Fig. 4 indeed confirms that the temperature at depth δ is a reasonable approximation.The temperature at z = −δ has been obtained by linear interpolation between levels, except when δ < 0.5∆z.
In the latter case, temperature T n 1 is selected.Note that, unlike α, β does change with temperature and does evolve during the integration.From Figs. 3 and 4, it is concluded that reasonable estimates can be made of α and β without actually solving the tri-diagonal matrix.Depth scale δ and the thickness of the top layer ∆z are crucial scales to characterize the temperature evolution of the top snow layer over a time step.

Simulations with the empirical formulation
With the empirical formulations for α and β, it is possible now to repeat the simulations of section 4. Instead of generating the fully implicit solution by solving the tri-diagonal matrix in the standard way, α and β are replaced by the empirical formulation between the elimination and back-substitution phase.If the formulation is perfect, the solution should be the same as the fully implicit solution.Results are shown in Fig. 5 for the skin temperature and the heat flux.Layer thicknesses of 0.2, 0.02 and 0.002 m are shown as different rows in Fig. 5.The figure confirms that the diurnal temperature cycle of the fully implicit solution (blue curve, IMPL) is well reproduced by the solution with parametrized α and β (black solid cure, IMPPAR).The differences between blue and black curves are very small.Finally, the scheme was further simplified by using the parametric form for α only and estimating β by putting it equal to 1 .The advantage is that no interpolation to z = −δ is needed, but that stability of the coupling is still maintained.However, it is clear that numerical errors are increased for thin snow layers (see dashed black curve).Such errors have to be seen in the context of other model errors, so the use of a parametrized α only, to ensure stability, may still be sufficient for many applications.

Discussion and conclusion
Numerical stability is a critical issue for atmospheric models that are coupled to a fast responding surface e.g. through a thin snow or ice layer.Very thin snow layers can occur in early winter after the first snow fall and during melt in spring.A fine discretization may also be desirable to allow for a fast response of the surface temperature to changes in radiation.Formal stability analysis confirms that unconditional stability can be achieved by a fully implicit coupling between atmosphere and surface.
Fully implicit coupling leads to a tri-diagonal problem in which atmosphere and surface are solved simultaneously.In practice, often so-called explicit flux coupling is applied: the atmospheric model uses the surface temperature of the previous time level to compute the surface heat flux, which is used later as boundary condition for the heat diffusion in the surface.
Explicit surface coupling puts stability limits on the thickness of the top snow layer and on the time step.Explicit flux coupling is often applied, because existing codes do not necessarily have sufficient modularity to support fully implicit coupling.
Although the atmosphere / surface heat diffusion leads to a single tri-diagonal matrix problem, one can also break it up in different steps.It is shown that the elimination part of the solver of the snow heat diffusion problem leads to a linear relation between surface temperature and surface heat flux.This relation can be used together with the atmosphere / surface interaction formulation to solve for the surface heat flux.
A simple method has been developed to approximate the coefficients in this linear relation.The coefficients are scaled with the characteristic scales of the diffusion equation.This makes the result universal and applicable to an arbitrary medium e.g.snow, ice or soil.The depth scale that characterizes the penetration of a perturbation over a time step, turns out to play a crucial role.In this paper the relevant empirical function is "measured" by solving the diffusion equation for a range of vertical resolutions and time steps.
Finally, the empirical functions are used to solve for the coupled diffusion problem and compared with the fully implicit computations.The results are very close.The advantage of the method is that the surface fluxes can be computed without calling any surface code, and behaves like explicit flux coupling.The only difference is that the surface heat flux expression has a damping term depending on the time step.This damping term is the result of the change of surface temperature related to the heat flux, and stabilizes the result.
The scaling argument used above only applies for a diffusion equation with constant properties of the medium.However, in reality there may be a profile of e.g.snow density as snow becomes more and more compact in deeper layers, or vertical resolution may be variable.The latter is numerically equivalent to a variable diffusion coefficient 1 .As a simple test, a case was selected where the profile of density is 150 kgm −3 at the surface, increases linearly to 250 kgm −3 at a depth of 0.5 m, and remains constant below 0.5 m.The characteristic depth is again computed as in section 5, and to non-dimensionalize, the snow properties are taken from the middle of the top snow layer.For this case the dimensionless α and characteristic temperature β are shown in Figs. 6 and 7.They are very close to the figures for constant snow properties (Figs. 3 and 4), which suggests that the sensitivity to snow properties is fairly small.In general, it is to be expected that the snow properties very close to the surface control the relation between flux and temperature over a short time step, because the penetration depth δ is small.We conclude that making an estimate of the relation between heat flux and surface temperature is a practical solution to support explicit flux coupling and to combine numerical stability for long time steps with a modular code structure.A formal stability analysis in Appendix B confirms unconditional stability of the proposed coupling method.The similarity framework makes the method applicable to any medium, e.g.snow, ice or soil.It is also worth noting that the method does not compromise conservation: the heat flux that is computed by the atmospheric model is later used by the surface model as boundary condition.
1 In fact the aerodynamic coupling between atmosphere and snow can be interpreted as a big jump in the properties of the medium Appendix A: Solving the tri-diagonal matrix equations The set of equations discussed in section 2 leads to the following tri-diagonal system where -Parametrized implicit flux coupling: the temperature at time level n + 1 is diagnosed as where α is defined in (21).Using ( 22) it can readily be seen that which shows that the parametrized implicit flux coupling can be interpreted as a limiter acting on the value of γ of the explicit flux coupling, indeed γ(1 + αλ t ) −1 ≤ γ.

B2 Matrix stability analysis
As shown in Appendix A, the Euler implicit scheme applied to the diffusion equation can be written in a general matrix form with where θ = 0 corresponds to the explicit flux coupling, θ = 1 to the implicit flux coupling, and θ = 1/2 to the parameterized implicit flux coupling.The general implicit scheme (B6) is stable if all the eigenvalues of the matrix M = A −1 B do not exceed 1 in magnitude.Therefore, the stability analysis requires the computation of the spectral radius of matrix M, i.e. its larger eigenvalue in magnitude.For γ ≥ 0 it can be shown that the smallest eigenvalue of A is larger or equal 2 to 1 meaning that this matrix is invertible for θ ∈ {0, 1/2, 1}.In Fig. 8, values of the spectral radius of M obtained over a range of values of γ 2 for the special cases θ = 0 and θ = 1/2, the eigenvalues λ A k of matrix A are given by λ A k = 1 + 2σ 1 + cos kπ N L for k = 1, . . ., N L; therefore and σ are shown for each coupling algorithm3 .Gray shaded areas coincide with regions where the spectral radius is larger than 1 thus indicating parameter values for which the corresponding scheme is unstable.From those results, the only algorithm that turns out to be conditionally stable is the explicit flux coupling whereas the implicit and parameterized implicit flux coupling are unconditionally stable.The results are thus consistent with the numerical experiments discussed in sections 4 and 6.
Empirically, it can be found that the stability condition for the explicit flux coupling roughly behaves like γ ≤ 2 + √ σ 1.1 (see figure 9).For the parametrized implicit flux coupling, γ is replaced by 3 ) 1/1.3 .As shown in figure 9, ∀σ ≥ 0, γ max ≤ 2 + √ σ 1.1 meaning that for the particular choice of f ( √ σ) given in ( 22) the parameterized implicit flux coupling is unconditionally stable because it always satisfies the stability constraint of the explicit flux coupling.

Figure 1 .Figure 2 .
Figure1.The numerical grid is defined by the position of the half levels, i.e. the thickness of the layers.The full levels are in the middle of the layers, i.e. zj = (z j−1/2 + z j+1/2 )/2.The surface is at z = 0.The bottom level is defined by the accumulated depth of all the layers.The temperature is defined on full levels and the heat fluxes are defined on half levels.

Figure 3 .
Figure 3. Dimensionless function f = α(K ρC/∆t) 1/2 as a function of x = δ/∆z.The circles and triangles are for different combinationsof ∆z and ∆t.The blue line is the asymptotic limit for small δ/∆z.The green curve is the empirical fit according to equation (22).

Figure 4 .Figure 5 .
Figure 4. Empirical estimates of parameter β as a function of the value found from the tri-diagonal solver.The red curve represents the estimate according to T n 1 and the blue curve is the temperature at z = −δ, also at the previous time level n.The symbols (connected by lines) indicate the successive time steps in the diurnal cycle.Results are plotted for vertical resolutions of 0.2, 0.02 and 0.002 m.

Figure 6 .Figure 9 .
Figure6.Dimensionless α as in Fig.3, but for non-uniform snow density.The snow density is 150 kgm −3 at the surface, increases linearly to 250 kgm −3 at a depth of 0.5 m, and remains constant below 0.5 m.

Table 1 .
List of parameters used in the idealized simulation of a snow layer