F. 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 "oi= aeA u*2i U"i |U"i |-1 Ei = aeA u*i Q*i (F:1) Hi = aeA CpA u*i *i L "i = - ffli oe Ti4 + ffLi L #; where the turbulent velocity scales are given by u*i = CD1=2i |U"i | Q*i = CEi |U"i | (qi) u*-1i (F:2) *i = CHi |U"i | (i) u*-1i ; where aeA is atmospheric surface density, CpA is the specific heat, oe = 5:67x10-8 W/m2 /K4 is the Stefan-Boltzmann constant, ffli is the emissivity of the interface and ffLi is the surface albedo for incident longwave radiation, L #. In (F.2) the differences U"i , qi 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, Zi, (Zi = ZA for i = 1, 2, 3), and to the stability, ii, appropriate to the particular interface, are : -2 CDi = 2 ln Zi___Zo- m i -1 Z -1 CEi = 2 ln Zi___Zo- m ln _i__e - s i Zi -1 Z -1 CHi = 2 ln Zi___Zo- m ln _i___h - s ; i Zi 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, ii. These functions as used in the coupler are: m (i) = s (i) = -5i i > 0 m (i) = 2ln[0:5(1 + X )] + ln[0:5(1 + X 2)] - 2tan-1 X + 0:5ss i < 0 s (i) = 2ln[0:5(1 + X 2)] i < 0 X = (1 - 16i)1=4 Above the atmospheric interfaces i = 1; 2 and 3 the stability parameter i * Q* j ii = _g_ZA______u*2i__+ ________i________-1 ; iv (Zv + qA ) F-1 where virtual potential temperature is computed as v = A (1 + Zv qA ), qA and A are the lowest level atmospheric humidity, and potential temperature, respectively, and Zv = (%(water)=%(air)) - 1 = 0:606. Since ii 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 ff (dir), and diffuse, ff (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, ffa (dif ). If direct and diffuse albedos, ffm (dir) and ffm (dif ) for m = 1; M different surfaces below an atmospheric grid cell are known, the multiple scattering coefficients, Rm = (1 - ffa (dif ) ffm (dif )) -1 are computed, and with f m the fractional coverage of each surface type, then the albedos are given by ffi ff(dif ) = (1 - F2 ) ( 1 - F2 =; ffa (dif )) ff(dir) = F1 ( 1 - ffa (dif ) ff(dif )) XM F1 = f m Rm ffm (dir) m=1 XM F2 = f m Rm (1 - ffm (dif )) m=1 Use of these albedos ensures that the solar radiation exiting the bottom of the at- mosphere is identical to the sum of M radiation calculations, each using an ffm (dif ) and ffm (dir). At present, however, the albedo calculations are greatly simplified by assuming ffa (dif ) = 0, such that Rm = 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 XM ff(dif ) = f m ffm (dif ) m=1 XM ff(dir) = f m ffm (dir) m=1 In general the partition of radiation among the M surfaces is a function of Rm , ffm (dif ), and ffm (dir), and hence of wavelength. Hence, correct partitioning would need to be performed for each wavelength band within the radiation transfer code of the atmo- spheric 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 ZA , and the ice surface temperature T1 = TC . The ice surface humidity is assumed to be saturated : q1 = ae-1A C5 exp (C6 =T1 ); where C5 = 640380: kg/m3 and C6 = -5107:4 K. The differences in (F.2) become U"1 = U"A - "U1 q1 = qA - q1 1 = A - T1 ; where the ice velocity, U"1 , is presently assumed to be negligible. Presently the coupler assumes constant and equal roughness lengths for ice : Zo1 = Zh1 = Ze1 = 0:04m: To start the iterative solver for the fluxes (F.1), i1 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 i1 , 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, ffl1 = 1, and the ice surface albedo for incident longwave radiation, ffL1 = 0:0. II. Fluxes Across the Atmosphere - Biosphere Interface (2) The fluxes across the Atmosphere - Biosphere Interface (2) are computed in the bio- sphere 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, q1 = 0:98 ae-1A C5 exp (C6 =T3 ); where the sea surface temperature is T3 = TH , and the factor 0.98 accounts for the salinity of the ocean. The differences (F.2) become U"3 = U"A - "U3 q3 = qA - q3 3 = A - T3 ; F-3 where the surface current, U"3 , is presently assumed to be negligible. The roughness length for momentum, Zo3 in meters, is a function of the atmospheric wind at 10 meters height, U10 : " -1 # Zo3 = 10 exp - C1____U+ C2 + C3 U10 ; 10 where C1 = 0:0027 m/s, C2 = 0:000142; and C3 = 0:0000764 m-1 s. The corresponding drag coefficient at 10m height and neutral stability is CN10 = C1 U1-10 + C2 + C3 U10 : The roughness length for heat, Zh3, is a function of stability, and for evaporation, Zo3, is a different constant: Zh3 = 2:2 x 10-9 m i3 > 0 = 4:9 x 10-5 m i3 0 Ze3 = 9:5 x 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, i3 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 i3 and U10 = UA . As with sea-ice, these coefficients are used to approximate the initial flux scales (F.2) and the first iteration begins with an updated i3 and calculations of m and s . The wind speed, UA , is then shifted to its equivalent neutral value at 10m height : i p ______CN Z j -1 U10 = UA 1 + _________ln10( _A___10- m (i3 )) : 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 i3 . 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, ffl3 = 1, and the water surface albedo for incident longwave radiation, ffL3 = 0:0. IV. Fluxes Across the Cryosphere - Hydrosphere Interface (4) Only the stress between the water and ice can be computed using bulk formulae sim- ilar to (F.1). The upper water column stability is assumed to be near neutral (i4 = 0:0; m (0) = 0:0), so that CD4 , and hence o"4 depend only on a constant roughness length and the velocity difference, Zo4 = 0:05m U"4 = U"C - U"H ; F-4 where U"C and UH" are, respectively, the ice velocity and surface current. The implied drag coefficient for Z4 = 10m is CD4 = 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 QH . The sense of this heat transfer is downwards, so QH 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 QH > 0. If there is no freezing, QH = aeH___CpH____1___(Tf__-_TH__)________t < 0 ; m is used to represent the potential of the upper ocean layer (thickness 1 ) to melt sea-ice, where tm 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 QH < 0 in any way. Instead the heat and freshwater fluxes associated with melting, both due to warm surface water QH , and atmospheric heating are computed within the cryosphere and passed to the Flux Coupler as H4 and F4 , respectively. The solar radiation passing through the bottom of the ice, S4 , is found in the cryosphere component by reducing the atmospheric solar, S1 , 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 Rn ; n = 1; N . The freshwater flux associated with runoff is passed to the ocean model as a vector F5 (m); m = 1; M , where each element of F5 is the runoff into one of the M ocean surface grid points. This flux is computed by a matrix multiplication in the coupler F5 = 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 Tm;n is the fraction of runoff from land point n that is transported to ocean grid point m. For an ocean grid of dimension Im longitudes and Jm latitudes, M = Im x Jm and m = i + (j - 1)Im at grid point (i, j). Similarly, N = In x Jn and n = i + (j - 1)In at grid point (i, j) of a land model with In longitudes and Jn 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 X Tm;n = 1:0 ; m=1 for all n. At present runoff is not transported by the coupler, i.e., Tn;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