The Role of Non‐Local Effects on Surface Sensible Heat Flux Under Different Types of Thermal Structures Over the Arctic Sea‐Ice Surface

The effects of atmospheric thermal structure on the surface energy flux are poorly understood over the Arctic sea‐ice surface. Here, we explore the mechanism of sensible heat exchange under different types of thermal structures over the Arctic sea‐ice surface by using data collected during the Multidisciplinary drifting Observatory for the Study of Arctic Climate expedition. The quadrant analysis indicates that strong surface temperature inversions below 100 m enhance non‐local effects on the positive (upward) sensible heat flux ( w′θ′‾ $\overline{{w}^{\prime }{\theta }^{\prime }}$ ) through entrainment of large eddies from the convective boundary layer aloft. However, strong surface inversions restrict the contributions of large eddies to the negative (downward) w′θ′‾ $\overline{{w}^{\prime }{\theta }^{\prime }}$ due to intensified surface stability. By inspecting the existing parameterization schemes, we found that the European Center for Medium‐Range Weather Forecasts Integrated Forecasting System scheme fails to predict the impacts of non‐local processes on the positive w′θ′‾ $\overline{{w}^{\prime }{\theta }^{\prime }}$ , and an adjustment term to correct the bias of parameterized w′θ′‾ $\overline{{w}^{\prime }{\theta }^{\prime }}$ is proposed.


Introduction
Arctic warming has attracted extensive attention, not only for its accompanying fast-decreasing sea ice but also for the more sensitive global climate system (Hartfield et al., 2018;Kwok, 2018;Landrum & Holland, 2020;Maure et al., 2023;Min et al., 2023).Processes occurring in the atmospheric boundary layer (ABL) play an essential role in understanding the ongoing climate changes in the Arctic (Wendisch et al., 2019;Wetzel & Bruemmer, 2011).The Arctic ABL is characterized by a predominant temperature inversion (Tjernström & Graversen, 2009).Some previous studies reported that the persistence of a shallow, single-layer stratocumulus at the top of the ABL was dominant over the Arctic Ocean (Leck et al., 1996;Tjernström et al., 2012;Uttal et al., 2002).However, Vüllers et al. (2021) found a high occurrence of multiple cloud layers with multi-layer temperature inversion structures in the central Arctic during the summer of 2018.Hence, the thermal structures of the Arctic ABL can be diverse.
The ABL and surface interact via the exchange of surface fluxes, and in the Arctic, this exchange is strongly influenced by the semipermanent sea ice coverage.To understand the evolution of the ABL, we need to clarify its interaction with the surface energy exchange (Persson, 2012).Many studies have confirmed that surface energy fluxes are not only generated by local processes but are also affected by non-local ones (De Bruin et al., 1999;Esau, 2004;Esters et al., 2021;Lohou et al., 2010;van de Boer et al., 2014).For the convective boundary layer (CBL), large-scale (∼10 3 m) coherent structures produce random gusts that crucially affect the surface fluxes (Grachev et al., 1998), such that the surface fluxes seem to depend on the global properties of the CBL (e.g., depth).Van de Boer et al. (2014) detected non-local effects on the surface humidity variance induced by entrainment processes and quantified entrainment influences on surface-layer observations.Meanwhile, nonlocal effects are also identified by experimental evidence in the stable boundary layer (SBL).For example, Zilitinkevich and Calanca (2000) related the free-flow static stability to the surface fluxes and proposed a nonlocal similarity based on summer measurements over the Greenland ice sheet.Zilitinkevich (2002) developed a theoretical model to present the non-local turbulent transport, accounting for the internal-wave-induced interaction between the long-lived nocturnal SBL and the free atmosphere.To our knowledge, it is widely accepted that non-local transport is affected by atmospheric thermal structures aloft, but its behavior under different types of thermal structures is complicated and poorly understood.As a result, using conventional localscale similarity theory (C.Liu et al., 2020) to calculate surface fluxes might be problematic until the non-local effects have been well considered (Mahrt & Bou-Zeid, 2020).
Scarce observations limit our understanding of the surface flux exchange over the Arctic sea-ice surface.Previous field campaigns focusing on atmosphere-ice-ocean interactions in the Arctic were mostly conducted during shortterm periods (Graham et al., 2017;Leck et al., 1996Leck et al., , 2001;;Tjernström et al., 2014).The last comprehensive campaign with an entire annual cycle was the Surface Heat Budget of the Arctic Ocean (SHEBA), conducted over 20 years ago (Uttal et al., 2002).Recently, to support the urgent need for understanding and modeling the rapidly changing Arctic atmosphere-ice-ocean system, an annual cycle expedition named the Multidisciplinary drifting Observatory for the Study of Arctic Climate (MOSAiC) was conducted from September 2019 to October 2020 in the central Arctic to collect a wealth of needed observational data (Shupe et al., 2022).The comprehensive and complementary atmospheric observational program during the expedition provides a unique opportunity to further reveal the mechanisms of surface energy exchange under different types of thermal structures.
The sensible heat flux quantifies one key component of the heat exchange between the sea-ice surface and atmosphere, and is an essential term of the sea-ice surface energy budget (Persson, 2012).In addition, it may also play an important role in interpreting climate changes in near-surface temperature over the Arctic sea-ice surface (Dai et al., 2019;Y. Liu et al., 2023;Screen & Simmonds, 2010).In this study, we aim to clarify the role of nonlocal effects on the surface sensible heat flux (w′θ′, here w′ and θ′ are turbulent fluctuations of vertical wind speed and potential temperature, respectively; the overbar denotes the time average, and the positive value indicates upward transport, and vice versa) under different types of thermal structures over the Arctic sea-ice surface and attempt to improve the parameterization of w′θ′ in atmospheric numerical models by taking non-local effects into account.The data and methods are described in Section 2.Then, the results and discussion are presented in Section 3. Finally, the summary and future perspectives are given in Section 4.

Observational Data
MOSAiC was based on and around the research icebreaker Polarstern as it drifted for a year with the sea ice across the central Arctic.Many details about the expedition are described in Shupe et al. (2022).In this study, w′θ′ and friction velocity (u * = w′u′ 2 + w′v′ 2 ) 1/ 4 , here u′ and v′ are turbulent fluctuations of horizontal velocity along the zonal and meridional directions, respectively) are obtained by the eddy-covariance (EC) technique using measurements from a u-Sonic-3 Cage MP anemometer (METEK GmbH, Germany) installed at a height of ∼6 m on a meteorological tower installed on the sea-ice surface (Figure S1 in Supporting Information S1) (Cox et al., 2023).The sampling frequency of the anemometer is 20 Hz, resampled to 10 Hz.The EC data processing done for this study includes error flag detection, despiking, true wind correction, coordinate rotation via double rotation, and block averaging over a 10-min period.The surface meteorology and flux data are collected during the period from 15 October 2019 to 18 September 2020.
In addition, a 6-hourly radiosonde program was undertaken by the Alfred Wegener Institute (AWI) and the US Department of Energy Atmospheric Radiation Measurement Program (Maturilli et al., 2021).The sounding data below ∼100 m altitude may be contaminated by the vessel.To avoid this contamination, a merged data product that combines the soundings with meteorological measurements from the tower installed on the sea ice is used to obtain information on the thermodynamic structures of the atmosphere (Dahlke et al., 2023).During the MOSAiC expedition, 1,484 radiosonde profiles in total were collected.In order to match the radiosonde measurements, the 10-min EC fluxes are averaged to hourly resolution within an hour after each radiosonde is released.
A hyperbolic hole (e.g., H = 1) is also applied to separate the relatively short-lived extreme events from the smallmagnitude fluctuations (i.e., |w′θ′| > H ⃒ ⃒ w′θ′ ⃒ ⃒ ).For each quadrant k and a hole size H, the average fluxes can be calculated as (Gao et al., 2018): where N is the number of samples for a 10-min period, and (2) When H = 0, the computed values of S for the four quadrants add up to 1 because the contributions at all scales combine to result in the total net flux.When H > 0, the computed values of S represent the fractions of flux contribution from events due to short-lived non-local motions at scales larger than H rather than the fluctuations at smaller scales.It is evident that the larger S with larger H corresponds to more contribution of large-scale eddies, and the surface fluxes are more likely to be affected by non-local processes.

Thermal Structure Characteristics
The prevalence of strong inversions is observed during the MOSAiC expedition, and the inversion strength is particularly intensified in winter (Figure 1a).To detect strong inversions for an individual potential temperature (θ) profile, the following procedures are carried out: (a) find the range of heights where the potential temperature lapse rate (∂θ/∂h) is greater than μ ∂θ/∂h + σ ∂θ/∂h (here μ and σ indicate the mean and standard deviation, respectively); (b) exclude inversion layers with a depth less than 30 m, which are likely caused by observational spikes; (c) if there are multiple strong inversions detected below the height of 2,000 m, only the inversion closest to the surface with the greatest impacts on the surface fluxes is retained.After these procedures, two types of θ profiles are recognized according to the location of the strong inversion.One is characterized by a strong inversion with a base below 100 m (Type1 or "surface inversion"); the other type has a strong inversion with a base height above 100 m (Type2).Type1 and Type2 profiles account for 33.7% and 66.3% of total samples, respectively.Figure S2 in Supporting Information S1 presents examples of these two types of θ profiles.
Figure 1b shows the statistics of θ-profile types in different seasons.The division of seasons is according to the near-surface atmospheric and sea-ice state as presented in Shupe et al. (2022), that is, defining the seasons of autumn, winter, spring, and summer roughly corresponds to the "freeze up" (from 11 October 2019 to 25 November 2019 and from 8 September 2020 to 1 October 2020), "winter" (from 26 November 2019 to 13 April 2020), "transition" (from 14 April 2020 to 25 May 2020), and "summer melt" (from 26 May 2020 to 7 September 2020) periods, respectively.It is found that the occurrence of strong surface inversions is relatively less frequent in each season.However, this condition becomes relatively more frequent in winter and summer.The increase of strong surface inversions during winter may be due to the radiative cooling of the sea-ice surface combined with relatively warm air advection aloft.The formation of strong surface inversions in summer can be attributed to the significant increase in air temperature coupled with the melting surface's constrained temperature (∼0°C) (Peng et al., 2023).
The distribution of w′θ′ for different types of θ profiles during the whole MOSAiC expedition is presented in Figure 1c.While negative w′θ′ is dominant under Type1 profiles, there are still some weak positive w′θ′ occurring under these strong surface inversions, indicating thin convective layers beneath the strong surface inversions.For Type2 cases, most samples (77.7%) are distributed in the range of 0.01 to +0.01 K m s 1 , and the probability of the positive w′θ′ is somewhat greater than that of the negative w′θ′ overall.In the following context, we will separately analyze the role of strong inversions in effecting the turbulent structure of the positive and negative w′θ′.

Turbulent Structures of Surface Sensible Heat Flux
It should be noted that non-local processes generally include advection, entrainment, or additional sources and sinks, however advection is assumed to be negligible in our study because the sea-ice surface is horizontally homogeneous over a large scale.Hence, the non-local processes in this study are regarded as entrainment for the CBL (De Bruin et al., 1993, 1999;Lohou et al., 2010;van de Boer et al., 2014) or internal-wave-induced interactions for the SBL (Zilitinkevich & Esau, 2005, 2007;Zilitinkevich et al., 2007).To further investigate the impacts of non-local effects on w′θ′ under different θ-profile conditions, the quadrant analysis is applied.The medians of quadrant analysis resulting from 10-min w′θ′ for the two thermal structure types are shown in Figure 2. To aid in interpreting Figure 2, we first provide a qualitative explanation of the meaning of each aspect of the figure.The different slopes of S diminish toward 0 with an increase in H, indicating that non-local effects play different roles in the four quadrants.Figure 2a presents the results for the positive w′θ′.
For these upward sensible heat fluxes, it is clear that the contributions of ejections and sweeps are both larger over the full range of H for Type1 cases than Type2 cases, especially for ejections.This indicates that Type1 cases have larger non-local contributions by large eddies to the positive w′θ′ compared with Type2 θ profiles.Hence, the above results indicate that strong surface inversions facilitate the contribution of large eddies to the positive w′θ′.In Type1 cases, the significant contributions of large eddies to the positive w′θ′ may occur through entrainment.In a CBL, entrainment causes large eddies to transport warm air from the strong surface inversion layer aloft down toward the surface, which can be verified from the quadrant analysis result of S 4H shown in Figure 2a.These warm air parcels then make considerable contributions to the positive w′θ′.However, for Type2 cases, the upper air carried by large eddies is relatively cold, and when it reaches the surface, it makes a smaller contribution to the positive w′θ′.Figure S3 in Supporting Information S1 gives a schematic illustration showing the effects of large eddies on the positive w′θ′ in the presence and absence of a strong surface inversion.
For the negative w′θ′ (Figure 2b), the contributions of ejections and sweeps are larger for a wide range of H when there is not a strong surface inversion (Type 2) compared to when there is a strong surface inversion (Type 1).This result means that a strong surface inversion inhibits the contributions of large eddies to downward sensible heat flux.In Type1 cases, the surface stability is relatively strong with a mean z/L (where z is the observation height and L is the Obukhov length) of 4.0, corresponding to a strong SBL.On the other hand, for Type2 cases, the mean z/L is equal to 0.18, corresponding to a weak SBL.As illustrated by Lan et al. (2018), turbulent eddies are confined to a thin layer under a strong SBL, leading to the dominant contribution of local gradients to the negative w′θ′.Under weak SBL conditions, the enhanced downward mixing of heat can be attributed to the development of turbulent eddies with larger scales.Hence, the contributions of large-scale eddies to the surface flux are more easily stimulated under weak SBL conditions.

Effects of Non-Local Processes on Parameterizing Surface Sensible Heat Flux
The Monin-Obukhov Similarity Theory (MOST) is widely used for parameterizing surface turbulent fluxes in numerical models.However, it is a local theory and does not incorporate the possible effects of non-local transport caused by large-scale eddies on the atmospheric surface layer (ASL) (C.Liu et al., 2020).As illustrated in Section 3.2, the contributions of large-scale eddies to w′θ′ cannot be overlooked and vary under different conditions.In this section, we discuss the effects of non-local mixing on the parameterization of surface w′θ′ under different types of θ-profile conditions.

Negative w′θ′
In the classic MOST framework, the negative w′θ′ can be parameterized as: Here, ∂u/∂z and ∂θ/∂z are the gradient of wind speed and potential temperature from the observation height (i.e., 6m height) to the surface, respectively; k (=0.4) and k T (=0.42) are the von Kármán constant with respect to wind and temperature, respectively; C u1 (=2.1) and C θ1 (=3.2) are known empirical constants (Högström, 1996); θ * is the potential temperature scale; and g (=9.8 m s 2 ) is the gravitational constant.In this framework, ∂u/∂z and ∂θ/ ∂z are required as input, then w′θ′ can be derived by iteratively solving Equations 3∼6.However, to consider nonlocal effects, Zilitinkevich and Calanca (2000) proposed an extended similarity theory for the stably stratified ASL as follows (only revised equations related to the heat flux are given here): Here, C θ2 is a dimensionless coefficient; Fi the inverse Froude number (higher values of Fi indicate more significant contributions from non-local transport); N is the Brunt-Väisälä frequency; and h i is the ABL height.In this study, the determination of h i follows Peng et al. (2023).To analyze the dependence of the surface scaling on Fi, Equation 7 can be rewritten as follows: First, we find the mean (±standard deviation) value of N is 2.5 (±0.7) × 10 2 s 1 .This value is slightly higher than the typical value of 10 2 s 1 , but is in accordance with values reported at high latitudes, for example, 3.0 × 10 2 s 1 for Arctic summer (Zilitinkevich & Calanca, 2000) and 2.0 × 10 2 s 1 over an Antarctic ice shelf (King, 1990).Once N is determined, Fi is calculated to quantify the effects of non-local transport accounting for the internal-wave-induced interaction between the long-lived SBL and the free atmosphere.The mean (±standard deviation) Fi in Type1 and Type2 cases is 8.6 (±5.6) and 11.6 (±6.2), respectively.Hence, the effects of non-local transport on the negative w′θ′ are less important in Type1 cases compared to Type2 cases.It is no surprise that the higher Fi occurs in Type2 cases because strong inversions always occur at high altitudes in Type2 cases.
Additionally, the dependence of the surface scaling on Fi (Figure 3) shows that the value of C θ2 in the best fitting curve for Type1 and Type2 cases is 0.5 and 0.6, respectively.That is to say, the surface scaling is more susceptible to non-local effects in Type2 cases.This also supports the results of the quadrant analysis and shows the suppressive effect of strong surface inversions on the contributions of large eddies to the negative w′θ′.The Fi and s θ given in Figure 3 are calculated using the EC observations rather than the parameterized results.This general result suggests that the coefficient C θ2 should be different under different types of θ-profile conditions by using Equations 7∼9.Overall, the values of C θ2 suggested in this study are smaller than the value of 1.0 recommended by Zilitinkevich and Calanca (2000) but still within the uncertainty of their results.In addition to the work of Zilitinkevich and co-authors on identifying non-local effects on the surface fluxes (Zilitinkevich & Calanca, 2000;Zilitinkevich & Esau, 2007;Zilitinkevich et al., 2002), Grachev et al. (2015) and Li et al. (2016) also related the Ozmidov scale (=ε 1/2 /N 3/2 , where ε is the dissipation rate of turbulent kinetic energy) to surface turbulent transport.

Positive w′θ′
In the CBL, Grachev et al. (1998) pointed out that the surface fluxes were affected by large-scale eddies encompassing the entire ABL.Therefore, the Deardorff velocity scale (ω * ) should be included in the parameterization of surface fluxes.The flux parameterization scheme used in the European Center for Medium-Range Weather Forecasts Integrated Forecasting System (ECMWF-IFS) is one example that uses this scaling; the details can be found in Supporting Information S1.Here, we use MOSAiC observations to evaluate the ECMWF-IFS scheme.
As shown in Figure S4 in Supporting Information S1, although ω * is included in the ECMWF-IFS scheme, the parameterized w′θ′ still presents a significant underestimation for cases of any type of θ profile.In addition, as presented in Equations S1∼S7 in Supporting Information S1, the parameterized w′θ′ has the same sign as the temperature gradient (θ s θ a ).However, Figure S4 in Supporting Information S1 indicates that there is a considerable number of negative parameterized w′θ′ values when the observed w′θ′ is positive.Such issues can be referred to as counter-gradient transport, that is, (θ s θ a ) is negative when the observed surface sensible heat flux is upward.Some previous studies have suggested that such underestimation and counter-gradient heat flux may be related to non-local transport (Blay-Carreras et al., 2014;Holtslag & Moeng, 1991;Smedman et al., 2007).Thus, these results suggest that the ECMWF-IFS scheme may not effectively explain the impacts of non-local transport on the positive w′θ′.Actually, the additional ω * term in the ECMWF-IFS scheme has very limited effects on the parameterization of surface fluxes compared with the traditional MOST framework (not shown) since ω * is much smaller than mean wind speed in Equation S3 in Supporting Information S1.Furthermore, the quadrant analysis shown in Section 3.2 reveals more significant contributions of non-local transport to the positive w′θ′ in Type1 cases than in Type2 cases.However, the contributions of non-local transport to the positive w′θ′ also cannot be ignored in Type2 cases from the perspective of parameterized results, as shown in Figure S4b in Supporting Information S1.Zilitinkevich et al. (1998) used the dimensionless minimum friction velocity (u * /ω * ) to parameterize w′θ′, which inspires us to revise the flux parameterization scheme further.Here, we use the dimensionless parameter ω 2 * / (u * θ * ) to describe the difference between the parameterized w′θ′ by the traditional MOST framework (w′θ′ p ) and the observed w′θ′ (w′θ′ o ).It is interesting that an exponential relationship is found between w′θ′ o w′θ′ p and ω 2 * / (u * θ * ) , as shown in Figure 4.The best-fit curve in Figure 4 is: (w′θ′ o w′θ′ p )/w′θ′ f = exp( 0.097(x + 33.66)) + 0.0015, (11) Here, w′θ′ f = 1 K m s 1 , is introduced to eliminate the dimension.Once the parameterized flux by the ECMWF-IFS scheme is solved, one can easily obtain the corrected w′θ′ (w′θ′ o ) according to: w′θ′ o = w′θ′ p + w′θ′ f × exp( 0.097(x + 33.66)) + 0.0015. (13) The independent SHEBA observations are used to further verify this relation.As shown in Figure S5 in Supporting Information S1, the parameterized w′θ′ has been significantly improved by introducing Equation 13 for the SHEBA data set.Accordingly, the new scheme produces an increase in the positive sensible heat flux by several W m 2 for the annual mean, which may partially correct the underestimated sensible heat flux over the Arctic sea-ice surface by the traditional MOST (Brunke et al., 2006).

Summary and Future Prospective
Previous studies have shown that surface sensible heat flux is affected by local-scale processes and ABL-scale structures.With increasing climate change in the Arctic, thermal structures in the ABL have attracted increasing attention.However, the effects of ABL thermal structures on w′θ′ over seaice surfaces are not thoroughly understood.Based on data collected during the MOSAiC expedition, the sources of w′θ′ under different types of thermal structures are investigated, and parameterization schemes that take non-local effects into account are evaluated and improved.
First, two types of θ profiles are identified according to the vertical location of strong inversions.The statistics show that the occurrence of a strong surface inversion accounts for 33.7% of all samples.It is interesting that upward w′θ′ is sometimes observed (∼25% of samples) beneath strong surface inversions, indicating the existence of thin and weak convective layers beneath the strong surface inversions.
Second, the quadrant analysis is used to analyze the flux contributions under the different types of θ-profile conditions.It is found that the effects of strong surface inversions on the negative and positive w′θ′ are contrasting.For the positive w′θ′, strong surface inversions facilitate the contributions of large eddies to the surface flux, likely via entrainment processes.However, strong surface inversions restrict the contributions of large eddies to the negative w′θ′, which is related to the near-surface stability.Overall, non-local effects on w′θ′ are found to be important.
Finally, non-local effects on the parameterization of w′θ′ are discussed.For the negative w′θ′, the solution of Zilitinkevich and Calanca (2000) with the empirical C θ2 of 0.5 and 0.6 agrees well with the MOSAiC observations for Type1 and Type2 θ-profile conditions, respectively.On the other hand, the ECMWF-IFS scheme does not effectively predict the impacts of non-local processes on the positive w′θ′.To correct the parameterized w′θ′ in a CBL, a scaling term of w′θ′ o w′θ′ p )/ w′θ′ f = exp( 0.097(x + 33.66)) + 0.0015, where x = ω 2 * / (u * θ * ) , is proposed.
This study improves the understanding of sensible heat flux on sea-ice surfaces under different types of atmospheric thermal structures.In future work, it is also worth considering how other ABL structures, for example, the dynamic structure characterized by low-level jet, affect the surface fluxes.In addition, we realize that the relation established to correct the bias of parameterized positive w′θ′ by the ECMWF-IFS scheme is empirical, although we have used the independent SHEBA observations to verify its applicability.The lack of a physical basis for this correction may hinder the promotion of this empirical relation to other regions.Thus, we plan to build a more physically interpretable parameterization scheme in the future.

Data Availability Statement
The radiosonde data is available on the PANGAEA archive via Dahlke et al. (2023).The MOSAiC surface flux and other meteorological data are available via Cox et al. (2023).The SHEBA data can be found in Andreas et al. (2007).
the tth sample falls into the quadrant k and satisfies |w′θ′| > H ⃒ ⃒ w′θ′ ⃒ ⃒ , otherwise I k|H (t) = 0.The fractions of flux contribution and time occupation of each quadrant k are calculated as: S k|H = w′θ′ k|H w′θ′ .

Figure 1 .
Figure 1.(a) The vertical and temporal evolution of potential temperature lapse rate, (b) the occurrence of various types of θ profile in different seasons, and (c) the distribution of w′θ′ for various types of θ profile.

Figure 2 .
Figure 2. Quadrant analysis of (a) the positive and (b) the negative w′θ′ under different types of θ-profile conditions. ∂u

Figure 3 .
Figure 3.The dependence of s θ on Fi in (a) Type1 and (b) Type2 cases.