## Abstract

In this study, a variational approach was employed to compute surface sensible heat flux over the Arctic sea ice. Because the variational approach is able to take into account information from the Monin–Obukhov similarity theory (MOST) as well as the observed meteorological information, it is expected to improve the pure MOST-based approach in computation of sensible heat flux. Verifications using the direct eddy-correlation measurements over the Arctic sea ice during the Surface Heat Budget of the Arctic Ocean (SHEBA) experiment period of 1997/98 show that the variational method yields good agreement between the computed and the measured sensible heat fluxes. The variational method is also shown to be more accurate than the traditional MOST method in the computation of sensible heat flux over the Arctic sea ice.

## 1. Introduction

Accurate representation of the surface sensible heat flux over the Arctic sea ice is important to characterize the energy transfer between the atmosphere and its underlying surface, and to better understand the interaction occurring in this coupled system (Cao et al. 2002). Precisely modeling the sensible heat flux, however, is a very challenging issue because computation of the surface sensible heat flux in a numerical weather prediction (NWP) model and a general circulation model (GCM) is still based on the Monin–Obukhov similarity theory (MOST), which has usually performed poorly for a very stable boundary layer (SBL). SBL research has been carried out over five decades (e.g., Grachev et al. 2005; Grachev et al. 2007a,b; Grachev et al. 2008), but there is not a systematic theory available yet to view a unified physical picture in a wide regime of atmospheric stability, even with the aid of finding another similarity regime characterized by the gradient-based scaling in very stable conditions (Sorbjan 2006a,c) and efforts on large-eddy simulations (Sorbjan 2006b).

To reduce uncertainties of MOST-based flux algorithms, we have employed the variational method in computations of sensible heat flux by making full use of the observed meteorological information and the information provided by MOST (e.g., Cao and Ma 2005; Cao et al. 2006). Through minimizing the differences between the computed and the observed meteorological variables, the variational method can sufficiently take into account the observed meteorological and surface conditions and can adjust the computed flux toward the optimal one (e.g., Xu and Qiu 1997).

The variational approach has recently been successfully applied to the computation of sensible heat flux over a flat homogeneous surface (Xu and Qiu 1997), a heterogeneous forest canopy (Cao and Ma 2005), and a large northern lake (Cao et al. 2006). It is shown in our previous studies that the variational method improves substantially computations of sensible heat flux over various underlying surfaces compared with the conventional MOST-based methods.

To our knowledge, little work has been done to employ the variational approach for computation of sensible heat flux over the Arctic sea ice. In this study, we intend to extend the variational approach for the computation of sensible heat flux over the Arctic sea ice. Also, we intend to verify the variational method–computed sensible heat flux against direct eddy-correlation measurements over the Arctic sea ice (Andreas et al. 1999; Persson et al. 2002; Grachev et al. 2005) and compare the variational method–computed sensible heat flux with that computed by the conventional flux–gradient method.

The next section presents the conventional flux–gradient and the variational method for the computation of sensible heat flux. Section 3 describes the datasets used in this work. In section 4, the computational results of sensible heat fluxes are presented and verified against the eddy-correlation-measured fluxes. Conclusions are given in section 5.

## 2. Description of methods

### a. Flux–gradient method

The heat flux can be computed using the flux–gradient method based on the following coupled equations (e.g., Yaglom 1977):

where *u* is the wind speed at the reference height *z*. Here, Δ*T* and Δ*q* are the differences of temperature and specific humidity at two heights (*z* and *z*_{0} or *z*_{0t} or *z*_{0q}). The von Kármán constant is *κ* = 0.4 and *u*_{*} is the friction velocity. The flux temperature scale and the flux specific humidity scale are denoted as *T*_{*} and *q*_{*}. They are related to the sensible heat flux (*H*) and latent heat flux (*λE*), respectively, by

where *ρ* is the air density, *c _{p}* is the specific heat of air at constant pressure,

*λ*E is the latent heat of evaporation, and

*E*is the evaporation rate. In Eqs. (2) and (3),

*F*= −

_{h}*u*

_{*}

*T*

_{*}and

*F*= −

_{q}*u*

_{*}

*q*

_{*}. In Eqs. (1)–(3)

*L*(e.g., Garratt 1994) is the Obukhov length, which is defined as

where *T* is the average temperature at two heights and *g* is the gravitational acceleration. The virtual temperature *T _{v}*, is defined as

It can be shown that *L* formulated in Eq. (5a) is equivalent to the following definition (e.g., Brutsaert 1982; Busch 1973; Cao and Ma 2005; Cao et al. 2006):

In Eqs. (1)–(3) *ψ _{m}*,

*ψ*, and

_{h}*ψ*are integral forms of the departure of wind speed, temperature, and moisture from their neutral values. For unstable conditions (

_{q}*z/L*< 0), they can be empirically expressed as (Businger et al. 1971)

where

The constants (*γ _{m}*,

*γ*) are equal to 13.1 and 7.9, respectively (Yaglom 1977). For stable conditions (

_{h}*z/L*> 0), the formulas suggested by Holtslag and deBruin (1988) and Beljaars and Holtslag (1991) are used:

where the constants *a*, *b*, *c*, and *d* are equal to 0.7, 0.75, 5, and 0.35, respectively. Also, we have tried the profile function derived from the Surface Heat Budget of the Arctic Ocean (SHEBA) data (Grachev et al. 2007a) and found that there are no notable differences between the profile function of Grachev et al. (2007a) and the one used in this study in terms of computation of sensible heat flux.

The roughness length *z*_{0} varies with ice and meteorological conditions. Based on their observations, Guest and Davidson (1991) found that values of *z*_{0} change from 3 × 10^{−4} m for very smooth ice to 3 × 10^{−2} m for very rough ice. The observations conducted by Andreas et al. (2004) indicated that variations of *z*_{0} over an ice could be in a magnitude of 10^{−7} to 10^{−2} m, depending upon the friction velocity and geographic locations. For the SHEBA data, we use *z*_{0} = 0.0005 m in this study. Also, we use *z*_{ot} = *z*_{0q} = 0.0001 m, as suggested by Andreas et al. (2004).

Equations (1) and (2) and the two unknown variables *u*_{*} and *F _{h}* consist of a well-determined system, where

*u*

_{*}and

*F*can be solved by an iterative procedure (e.g., Cao and Ma 2005; Cao et al. 2006). Because we are interested in surface sensible heat flux, the observed wind speed and temperature difference at the lowest level (2.2 m) are used for computations of the flux–gradient method.

_{h}### b. Variational method

Unlike the flux–gradient method, the variational method is flexible in terms of fully using observed meteorological information over an underlying surface, such as an ice, through a procedure to minimize the difference between observed and computed meteorological variables. This difference is defined as a cost function *J*, which measures errors between the calculated and observed wind speed, temperature difference, and moisture difference (e.g., Cao and Ma 2005; Cao et al. 2006):

Different from Xu and Qiu (1997), we do not use any information of surface energy balance as a constraint in the cost function (6) simply because there is no soil heat flux available in SHEBA observations. In Eq. (6), *u _{i}*

^{obs}, Δ

*T*

_{i}^{obs}, and Δ

*q*

_{i}^{obs}are the observed wind speed, temperature difference, and moisture difference between

*z*(at levels of 2.2, 3.2, 5.1, and 8.9 m) and

*z*

_{ot}or

*z*

_{0q}, respectively. In the cost function

*u*, Δ

*T*, and Δ

*q*are computed using the flux–gradient relations (1)–(3). Here,

*W*,

_{ui}*W*, and

_{Ti}*W*are weights for the wind speed, the temperature difference, and the moisture difference, respectively.

_{qi}The weights are, in principle, chosen to be inversely proportional to their respective maximally tolerated observation error variances (Cao and Ma 2005; Cao et al. 2006). Because of uncertainties in measurement errors, the weights are usually difficult to determine precisely, and they can be specified empirically (Daley 1996). In this study, the estimated values of weights are set as follows: *W*_{u1} = 12 (m^{−2} s^{2}), *W*_{T1} = 20 (K^{−2}), and *W*_{q1} = 10 (kg kg^{−1}). By choosing *W _{ui}* =

*W*

_{u1}/5,

*W*=

_{Ti}*W*

_{T1}/5, and

*W*=

_{qi}*W*

_{q1}/5 (

*i*= 2, 3, 4),

*u*

_{i}^{obs}, Δ

*T*

_{i}^{obs}, and Δ

*q*

_{i}^{obs}are used as a supplemental data source (e.g., Xu and Qiu 1997). These weights are within a range of the maximally tolerated observation error variances, which lead to better results for sensible heat flux computation. As mentioned in section 4, the variational computations are not very sensitive to the choice of the weights.

The cost function *J*, which is considered as the function of (*u*_{*}, *F _{h}*) in the variational method, is minimized to obtain optimal estimates of

*u*

_{*}and

*F*. This necessitates that the gradients of

_{h}*J*with respect to unknown variables

*u*

_{*}and

*F*become zero:

_{h}The analytic formulations of these gradient components are given in the appendix.

A quasi-Newton method is employed to find the minimum of the cost function *J*. This method needs an iterative procedure to compute the cost function (6) and its gradients (7), that is, (A1)–(A8), and to determine the minimum of *J* along a search direction. For Eq. (6), the inputs are the measured *u*, Δ*T*, Δ*q*, and the outputs are calculated *u*, Δ*T*, and Δ*q*, and estimated *u*_{*} and *F _{h}*. The following iterative procedure is performed to compute the cost function and its gradient and to derive

*u*

_{*}and

*F*:

_{h}Set up initial guesses of unknowns, say,

*u*_{*}= 0.3 m s^{−1}and*F*= 0.03 K m s_{h}^{−1}.Calculate

*L*from Eq. (5c).Calculate the cost function (6) and its gradients with respect to

*u*_{*}from Eqs. (A1) and (A3)–(A5) and with respect to*F*from Eqs. (A2) and (A6)–(A8)._{h}Perform the quasi-Newton method to search for zeros of the gradients of

*J*, so as to minimize the cost function*J*; the expected*u*_{*}and*F*are reached at the minimum of the cost function._{h}Repeat steps 2–5 until the procedure converges.

The convergence criterion is set to 10^{−4}, and the convergence requirement on searching zeros of gradients is set to 10^{−7}. Note that Δ*q* can be computed when *u*_{*} is obtained through the variational procedure; *L* and *F _{q}* are obtained through the iterative procedure.

It should be mentioned that the flux observations are not used in the variational procedure, which makes the method robust. The measured fluxes are, however, used for verification purposes to evaluate how good the variational method is in comparison with other methods such as the flux–gradient method. This requires the observed flux not being used in the flux–gradient method as well.

In the variational approach, Eqs. (1)–(3) and the unknown variables *u*_{*} and *F _{h}* form an overdetermined system, where

*u*

_{*}and

*F*are solved through the procedure to minimize the difference between observed and computed meteorological variables. The overdetermined problem is suitable for the variational method but not for the flux–gradient method (Cao et al. 2006). For the overdetermined system with two unknowns (

_{h}*u*

_{*}and

*F*) and three equations (related to u, Δ

_{h}*T*, and Δ

*q*), the variational method can improve the flux–gradient method via making use of additional information of Δ

*q*. Furthermore, the variational method is different from the flux–gradient approach in terms of algorithms to minimize the errors.

## 3. Data

The dataset used in this study was collected from October 1997 to October 1998 over the Arctic sea ice during the SHEBA observation period (Andreas et al. 1999; Persson et al. 2002). The turbulent and meteorological variables were measured at five nominal heights of 2.2, 3.2, 5.1, 8.9, and 18.2 m. The observed variables and measurement method are described in Andreas et al. (1999) and Persson et al. (2002). The data collected in the field campaign include sensible heat flux, latent heat flux, wind speed, temperature difference, and moisture difference. The sensible heat fluxes were measured using the eddy-correlation system. Since we are interested in surface sensible heat flux, the observed sensible heat flux at the lowest level (2.2 m) is used for a verification purpose for both methods. The wind speed, temperature difference, and moisture difference observed at four low levels are used to retrieve sensible heat flux, whereas variables observed at the upper level (18.2 m) are not used as a result of icing and other problems (Persson et al. 2002). All the data used in this study have a 1-h time interval. Readers are referred to Andreas et al. (1999) and Persson et al. (2002) for details.

After data quality control, the size of the observational data collected over the Arctic sea ice was reduced because of the data missing and the requirement of simultaneous availability of all data involved in computations. However, there are still about 2700 data points available, which are quite representative of meteorological conditions during the SHEBA observation period. In this study, 1 January 1997 is defined as Julian day 1.

## 4. Results

The flux–gradient and the variational methods described in section 2 are employed for computing sensible heat fluxes, which are verified against the direct eddy-correlation measured sensible heat fluxes.

The correlation diagrams between the observed and the computed sensible heat flux are shown in Fig. 1. The agreements between the measured and the variational method–computed sensible heat flux are good (Fig. 1). The correlation coefficient (*r*_{1}) between the observed and the computed sensible heat flux is 0.83 (Table 1). For the flux–gradient method, the correlation coefficient (*r*_{2}) reaches 0.75. Comparisons between Fig. 1a and Fig. 1b indicate that the variational method improves the flux–gradient method in computation of sensible heat flux, although the flux–gradient method captures the average tendency better (see trend lines in Fig. 1) as a result of the frequent occurrence of near-neutral conditions during the SHEBA observations. The improvement of the variational method over the flux–gradient approach is mainly achieved through using observational data at multiple levels. This overdetermined system is allowed for the variational method but not for the flux–gradient approach (Cao et al. 2006).

To understand if the difference between two correlation coefficients is statistically significant, we compute a statistic *z* (Dowdy and Wearden 1983),

where *z*_{r1} and *z*_{r2} are correlation coefficients, with a near-normal distribution, transferred from original correlation coefficients *r*_{1} and *r*_{2}. Here *n*_{2} and *n*_{1} (*n*_{2} = *n*_{1} = 2676) are the sampling size associated with transferred correlation coefficients *z*_{r1} and *z*_{r2}. Under the null hypothesis *H _{o}*, Eq. (8) has a normal distribution with

*n*

_{2}+

*n*

_{1}− 2 degrees of freedom. Our computation shows that the statistic

*z*is equal to 7.5 (>

*z*

_{0.005/2}= 2.9) with the statistically significant level of 99.5%, indicating that the correlation coefficient for the variational method is significantly better than the correlation coefficient for the flux–gradient method.

The further quantitative evaluation on the performance of two methods can be realized by computing the mean absolute error (MAE), defined as

where *F _{ic}* and

*F*are

_{io}*i*th (

*i*= 1, … ,

*N*, where

*N*is the number of the total time level of measurements) calculated and observed sensible heat fluxes. As demonstrated in Table 1, the flux–gradient method has an MAE of 4.5 W m

^{−2}, whereas the variational method has an MAE of 4.0 W m

^{−2}.

The statistical difference between two MAE can be determined using a *t* test through computation of statistic *t* (Dowdy and Wearden 1983):

where *n*_{2} and *n*_{1} (*n*_{2} = *n*_{1} = 2676) indicate the size of the *y*_{2} and *y*_{1}, *y*_{2} and *y*_{1} are the sample MAE for the flux–gradient and the variational methods, respectively, and *s*_{1} and *s*_{1} are the estimated standard deviations. Under the null hypothesis *H _{o}*, Eq. (10) has a

*t*distribution with

*n*

_{2}+

*n*

_{1}– 2 degrees of freedom. Based on the calculation, the statistic

*t*is equal to 4.3 (>

*t*

_{0.005/2}= 2.6) with the statistically significant level of 99.5%, indicating that the MAE of the variational method is significantly better than that of the flux–gradient method.

The variational approach also captures observed diurnal variations of the sensible heat flux well, especially for stable conditions. The example presented in Fig. 2 is representative of typical diurnal variations of the sensible heat flux over the Arctic sea ice. It should be mentioned that since some of the data were eliminated through the quality control procedure, the result shown in Fig. 2 is the sample that passed the quality control. Figure 2 illustrates a case that the observed sensible heat flux mainly decreased with time in the early morning and increased afterward on Julian day 315. The variational method–calculated sensible heat flux (MAE = 1.0 W m^{−2}) better resembles the diurnal variation and the magnitude of the observed sensible heat flux than the flux–gradient method (MAE = 1.6 W m^{−2}), especially during the early morning (0–4 h) and the afternoon (15–20 h). Over the period of 15–20 h, the flux–gradient method calculated flux even shows opposite signs to the observed sensible heat flux.

Sensitivity experiments have been performed to examine the effects of the weights in the cost function on the variational method–computed sensible heat flux. Similar to Cao and Ma (2005) and Cao et al. (2006), the results indicate that the variational computations are not very sensitive to the choice of the weights (not shown).

It should be pointed out that the variational method needs more iterations to reach a convergence for a given data point than the flux–gradient method.

## 5. Conclusions

The variational method is applied to computations of the sensible heat flux over the Arctic sea ice in this work. Since the variational method is capable of fully taking into account information of the existing MOST and the measured meteorological conditions over the underlying surface, it improves the conventional MOST-based flux–gradient method in computations of sensible heat flux over the Arctic sea ice.

It is shown that the variational method yields higher correlations between the computed and the measured sensible heat flux than the flux–gradient method, and that the variational method is more accurate than the conventional method. These improvements are mainly a result of more observed information over the Arctic sea ice being used by the variational method than the flux–gradient method.

It is demonstrated in our previous and present studies that the variational method is very useful in computations of sensible heat flux over various underlying surfaces (e.g., forest, water, and ice) and over both mid- and high latitudes. The variational method developed in our studies is planned to be applied for computation of sensible heat flux in numerical weather prediction models and general circulation models, using past and/or current weather and climate observations.

## Acknowledgments

We thank the Atmospheric Surface Flux Group (ASFG) members for their efforts in collecting and processing the data used in this work and three anonymous reviewers for their valuable comments and suggestions.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

### APPENDIX

#### Formulations for the Gradient Components

The detailed analytic expression of the gradient components in Eq. (7) can be derived as follows. From Eq. (6), we have

The derivative of the wind speed with respect to *u*_{*} is deduced from Eq. (1):

where

The derivatives of *ψ _{m}* with respect to

*u*

_{*}at

*z*

_{0}can be derived in the same manner.

Based on Eq. (2), the derivative of Δ*T* with respect to *u*_{*} can be expressed as

where

The derivative of *ψ _{h}* at

*z*

_{0t}with respect to

*u*

_{*}can be derived in the same manner.

The derivative of the specific humidity with respect to *u*_{*} can be obtained from Eq. (3):

where

The derivative of *ψ _{q}* at

*z*

_{0q}with respect to

*u*

_{*}can be derived in the same manner. Here,

*F*can be obtained by an iterative procedure.

_{q}The derivatives of *u*, Δ*T*, and Δ*q* with respect to the heat flux *F _{h}* are given by

where

where

and

where

The derivatives of *ψ _{m}*,

*ψ*, and

_{h}*ψ*with respect to

_{q}*F*at

_{h}*z*

_{0},

*z*

_{0t}, and

*z*

_{0q}can be derived in the same manner.

## Footnotes

*Corresponding author address:* Dr. Zuohao Cao, Meteorological Service of Canada, 4905 Dufferin Street, Toronto, ON M3H 5T4, Canada. Email: zuohao.cao@ec.gc.ca