F. \underline{Flux Formulae}

Flux computations in the Flux Coupler are denoted by a circle in Fig. 1. Many of the fluxes across the interfaces i = 1, 2, 3 and 4 are calculated from bulk formulae and general expressions are

\begin{aligned}τ⃗_{i} & = ρ_{A} u_{i}^{∗2} ΔU⃗_{i} |ΔU⃗_{i}|^{−1} \\ E_{i} & = ρ_{A} u_{i}^{∗} Q_{i}^{∗} \\ H_{i} & = ρ_{A} Cp_{A} u_{i}^{∗} θ_{i}^{∗} \\ L ↑_{i} & = − ϵ_{i} σ T_{i}^{4} + α_{i}^{L} L ↓, \end{aligned}  (F.1)

where the turbulent velocity scales are given by

\begin{aligned}u_{i}^{∗} & = CD_{i}^{1/2} |ΔU⃗_{i}| \\ Q_{i}^{∗} & = CE_{i} |ΔU⃗_{i}| (Δq_{i}) u_{i}^{∗−1} \\ θ_{i}^{∗} & = CH_{i} |ΔU⃗_{i}| (Δθ_{i}) u_{i}^{∗−1}, \end{aligned}  (F.2)

where ρ_{A} is atmospheric surface density, Cp_{A} is the specific heat, σ = 5.67×10^{−8}W/m^{2}/K^{4} is the Stefan-Boltzmann constant, ϵ_{i} is the emissivity of the interface and α_{i}^{L} is the surface albedo for incident longwave radiation, L ↓. In (F.2) the differences ΔU⃗_{i}, Δq_{i} and Δθ_{i} are defined at each interface in accord with the convention of fluxes being positive down.

The transfer coefficients in (F.2), shifted to a height, Z_{i}, (Z_{i} = Z_{A} for i = 1, 2, 3), and to the stability, ζ_{i}, appropriate to the particular interface, are :

\begin{aligned}CD_{i} & = κ^{2} [ln (\frac{Z_{i}}{Z_{i}^{o}} ) − ψ_{m}]^{−2} \\ CE_{i} & = κ^{2} [ln (\frac{Z_{i}}{Z_{i}^{o}} ) − ψ_{m}]^{−1} [ln (\frac{Z_{i}}{Z_{i}^{e}} ) − ψ_{s}]^{−1} \\ CH_{i} & = κ^{2} [ln (\frac{Z_{i}}{Z_{i}^{o}} ) − ψ_{m}]^{−1} [ln (\frac{Z_{i}}{Z_{i}^{h}} ) − ψ_{s}]^{−1} , \end{aligned}

where κ = 0.4 is von Karman’s constant and the integrated flux profiles, ψ_{m} for momentum and ψ_{s} for scalars, are functions of the stability parameter, ζ_{i}. These functions as used in the coupler are:

\begin{aligned}ψ_{m}(ζ) & = ψ_{s}(ζ) = −5ζ  ζ > 0 \\ ψ_{m}(ζ) & = 2ln[0.5(1 + X)] + ln[0.5(1 + X^{2})] − 2tan^{−1}X + 0.5π ζ < 0 \\ ψ_{s}(ζ) & = 2ln[0.5(1 + X^{2})]  ζ < 0 \\ X & = (1 − 16ζ)^{1/4} \end{aligned}

Above the atmospheric interfaces i = 1, 2 and 3 the stability parameter

ζ_{i} = \frac{κ g Z_{A}}{u_{i}^{∗2}} (\frac{θ_{i}^{∗}}{θ_{v}} + \frac{Q_{i}^{∗}}{(Z_{v}^{−1} + q_{A})}) ,

F-1 where virtual potential temperature is computed as θ_{v} = θ_{A}(1 + Z_{v}q_{A}), q_{A} and θ_{A} are the lowest level atmospheric humidity, and potential temperature, respectively, and Z_{v} = (ϱ(water)/ϱ(air)) − 1 = 0.606. Since ζ_{i} is itself a function of the turbulent scales (F.2), and hence the fluxes, an iterative procedure is generally required to solve (F.1).

In addition to surface fluxes, the atmospheric model requires effective surface albedos for both direct α (dir), and diffuse, α (dif), radiation at each wavelength. They are used in a single call to the computationally demanding atmospheric radiation routines. This call gives downward atmospheric albedo for diffuse radiation, α_{a}(dif). If direct and diffuse albedos, α^{m}(dir) and α^{m}(dif) for m = 1, M different surfaces below an atmospheric grid cell are known, the multiple scattering coefficients, R^{m} = (1 − α_{a}(dif) α^{m}(dif))^{−1} are computed, and with f^{m} the fractional coverage of each surface type, then the albedos are given by

\begin{aligned}α(dif) & = (1 − F_{2}) / (1 − F_{2}/, α_{a}(dif)) \\ α(dir) & = F_{1} (1 − α_{a}(dif) α(dif)) \end{aligned}

F_{1} = ∑_{m=1}^{M} f^{m} R^{m} α^{m}(dir)

F_{2} = ∑_{m=1}^{M} f^{m} R^{m} (1 − α^{m}(dif))

Use of these albedos ensures that the solar radiation exiting the bottom of the atmosphere is identical to the sum of M radiation calculations, each using an α^{m}(dif) and α^{m}(dir). At present, however, the albedo calculations are greatly simplified by assuming α_{a}(dif) = 0, such that R^{m} = 1. This albedo then does not need to be passed from the atmosphere to the coupler, and the coupler simply computes the effective albedos to be passed to the atmosphere as

\begin{aligned}α(dif) & = ∑_{m=1}^{M} f^{m} α^{m}(dif) \\ α(dir) & = ∑_{m=1}^{M} f^{m} α^{m}(dir) \end{aligned}

In general the partition of radiation among the M surfaces is a function of R^{m}, α^{m}(dif), and α^{m}(dir), and hence of wavelength. Hence, correct partitioning would need to be performed for each wavelength band within the radiation transfer code of the atmospheric model. This procedure is greatly simplified by partitioning the net solar radiation within the Flux Coupler.

F-2 I. Fluxes Across the Atmosphere - Cryosphere Interface (1)

The state variables required to calculate the fluxes (F.1) for i = 1, are the atmospheric wind velocity, potential temperature, density and specific humidity in (kg/kg) at a height Z_{A}, and the ice surface temperature T_{1} = T_{C}. The ice surface humidity is assumed to be saturated :

q_{1} = ρ_{A}^{−1} C_{5} exp(C_{6}/T_{1}),

where C_{5} = 640380. kg/m^{3} and C_{6} = −5107.4 K. The differences in (F.2) become

\begin{aligned}ΔU⃗_{1} & = U⃗_{A} − U⃗_{1} \\ Δq_{1} & = q_{A} − q_{1} \\ Δθ_{1} & = θ_{A} − T_{1} , \end{aligned}

where the ice velocity, U⃗_{1}, is presently assumed to be negligible.

Presently the coupler assumes constant and equal roughness lengths for ice :

Z_{1}^{o} = Z_{1}^{h} = Z_{1}^{e} = 0.04m.

To start the iterative solver for the fluxes (F.1), ζ_{1} is set to zero : ψ_{m} = ψ_{s} = 0, and the neutral coefficients are used to approximate the flux scales (F.2). The iteration begins by computing ζ_{1}, finding the non-neutral transfer coefficients and updating the flux scales. After two such iterations there is sufficient convergence to compute the fluxes (F.1). Also in (F.1), the reflected downward incident longwave radiation is simply accounted for by assuming an emissivity, ϵ_{1} = 1, and the ice surface albedo for incident longwave radiation, α_{1}^{L} = 0.0.

II. Fluxes Across the Atmosphere - Biosphere Interface (2)

The fluxes across the Atmosphere - Biosphere Interface (2) are computed in the biosphere model. Their formulae can be written in the forms of (F.1) and are detailed in Bonan (1996), “A Land Surface Model (LSM version 1.0) For Ecological, Hydrological, and Atmospheric Studies : Technical Description and User’s Guide”.

III. Fluxes Across the Atmosphere - Hydrosphere Interface (3)

At the atmosphere-hydrosphere interface the near-surface air is also assumed to be saturated,

q_{1} = 0.98 ρ_{A}^{−1} C_{5} exp(C_{6}/T_{3}),

where the sea surface temperature is T_{3} = T_{H}, and the factor 0.98 accounts for the salinity of the ocean. The differences (F.2) become

\begin{aligned}ΔU⃗_{3} & = U⃗_{A} − U⃗_{3} \\ Δq_{3} & = q_{A} − q_{3} \\ Δθ_{3} & = θ_{A} − T_{3}, \end{aligned}

F-3 where the surface current, U⃗_{3}, is presently assumed to be negligible.

The roughness length for momentum, Z_{3}^{o} in meters, is a function of the atmospheric wind at 10 meters height, U_{10}:

Z_{3}^{o} = 10 exp − [κ (\frac{C_{1}}{U_{10}} + C_{2} + C_{3}U_{10})^{−1}] ,

where C_{1} = 0.0027 m/s, C_{2} = 0.000142, and C_{3} = 0.0000764 m^{−1}s. The corresponding drag coefficient at 10m height and neutral stability is

C_{10}^{N} = C_{1} U_{10}^{−1} + C_{2} + C_{3} U_{10} .

The roughness length for heat, Z_{3}^{h}, is a function of stability, and for evaporation, Z_{3}^{o}, is a different constant:

Z_{3}^{h} = 2.2 × 10^{−9}m ζ_{3} > 0
= 4.9 × 10^{−5}m ζ_{3} ≤ 0
Z_{3}^{e} = 9.5 × 10^{−5}m .

Since the roughness lengths are neither constant nor equal, the iterative solution for the fluxes (F.1) differs somewhat than that for sea–ice. First, ζ_{3} is set incrementally greater than zero when the air–sea temperature difference suggests stable stratification, otherwise it is set to zero. In either case, ψ_{m} = ψ_{s} = 0, and the initial transfer coefficients are then found from the roughness lengths at this ζ_{3} and U_{10} = U_{A}. As with sea–ice, these coefficients are used to approximate the initial flux scales (F.2) and the first iteration begins with an updated ζ_{3} and calculations of ψ_{m} and ψ_{s}. The wind speed, U_{A}, is then shifted to its equivalent neutral value at 10m height :

U_{10} = U_{A} (1 + \frac{\sqrt{C_{10}^{N}}}{κ} ln(\frac{Z_{A}}{10} − ψ_{m}(ζ_{3})) )^{−1} .

This wind speed is used to update the transfer coefficients and hence the flux scales. The second and final iteration begins with another update of ζ_{3}. The final flux scales then give the fluxes calculated by (F.1). As for the sea–ice interface (1), the reflected downward incident longwave radiation is simply accounted for by assuming an emissivity, ϵ_{3} = 1, and the water surface albedo for incident longwave radiation, α_{3}^{L} = 0.0.

IV. Fluxes Across the Cryosphere - Hydrosphere Interface (4)

Only the stress between the water and ice can be computed using bulk formulae similar to (F.1). The upper water column stability is assumed to be near neutral (ζ_{4} = 0.0; ψ_{m}(0) = 0.0), so that CD_{4}, and hence τ⃗_{4} depend only on a constant roughness length and the velocity difference,

\begin{aligned}Z_{4}^{o} & = 0.05m \\ ΔU⃗_{4} & = U⃗_{C} − U⃗_{H} , \end{aligned}

F-4 where U⃗_{C} and U⃗_{H} are, respectively, the ice velocity and surface current. The implied drag coefficient for Z_{4} = 10m is CD_{4} = 0.0057.

However, with the present sea–ice model, the stress is not computed in this fashion, but is given by a linear function of ΔU⃗_{4}, as described in Section G.

The heat and freshwater fluxes between the cryosphere and hydrosphere due to the freezing and melting of sea-ice are computed, respectively, in the hydrosphere (see ocean model documentation) and cryosphere (see ice model documentation) component models. Each column of the hydrospheric component needs to be scanned after each time step to see if any ice is formed. The cumulative heat of fusion associated with the freezing is output as Q_{H}. The sense of this heat transfer is downwards, so Q_{H} is positive. The heat and freshwater fluxes associated with freezing are applied directly by the hydrospheric model and not passed to the coupler other than by the cumulative Q_{H} > 0.

If there is no freezing,

Q_{H} = \frac{ρ_{H} Cp_{H} Δ_{1} (T_{f} − T_{H})}{t_{m}} < 0 ,

is used to represent the potential of the upper ocean layer (thickness Δ_{1}) to melt sea–ice, where t_{m} is a melting time scale. Since the hydrosphere itself does not know how much ice is available to be melted it does not respond to Q_{H} < 0 in any way. Instead the heat and freshwater fluxes associated with melting, both due to warm surface water Q_{H}, and atmospheric heating are computed within the cryosphere and passed to the Flux Coupler as H_{4} and F_{4}, respectively.

The solar radiation passing through the bottom of the ice, S_{4}, is found in the cryosphere component by reducing the atmospheric solar, S_{1}, to account for heat absorption by the ice and snow.

V. Fluxes Across the Biosphere - Hydrosphere Interface (5)

The only important direct exchange between the biosphere and hydrosphere (across interface 5) is river runoff. The runoff from each of the N land model grid points is passed to the coupler as a vector R_{n}, n = 1, N. The freshwater flux associated with runoff is passed to the ocean model as a vector F_{5}(m), m = 1, M, where each element of F_{5} is the runoff into one of the M ocean surface grid points. This flux is computed by a matrix multiplication in the coupler

F_{5} = T R ,

where T is an (M by N) runoff transport matrix that is specified in the coupler initialization and remains constant in time. Matrix element T_{m,n} is the fraction of runoff from land point n that is transported to ocean grid point m.

For an ocean grid of dimension I_{m} longitudes and J_{m} latitudes, M = I_{m} × J_{m} and m = i + (j − 1)I_{m} at grid point (i, j). Similarly, N = I_{n} × J_{n} and n = i + (j − 1)I_{n} at grid point (i, j) of a land model with I_{n} longitudes and J_{n} latitudes. Note, runoff from a

F-5 land gridpoint can be transported to different ocean basins, including marginal seas where it can be made to balance the local evaporation and precipitation. River runoff can be distributed over a number of ocean grid points. Conservation of freshwater is satisfied by the condition

∑_{m=1}^{M} T_{m,n} = 1.0 ,

for all n.

At present runoff is not transported by the coupler, i.e., T_{n,m} = 0 for all n and m. The freshwater contained in the active ocean is conserved in the coupler by multiplying the precipitation by the uniform factor needed to balance the global average evaporation every coupling period, see section G.

F-6