Nonlinear Processes in Geophysics Scale-dependency of surface fluxes in an atmospheric mesoscale model : effect of spatial heterogeneity in atmospheric conditions

We examined the nonlinear effect of spatial heterogeneity in atmospheric conditions on the simulation of surface fluxes in the mesoscale model, MM5 by testing their scale-invariance from a tower footprint to regional scales. The test domain was a homogeneous shortgrass prairie in the central part of the Tibetan Plateau with an eddy-covariance flux tower at the center. We found that the spatial variability resulting from changing distribution of clouds and precipitation in the model domain affected radiative forcing at the ground surface, thereby altering the partitioning of surface fluxes. Consequently, due to increasing spatial variability in atmospheric conditions, the results of MM5 did not produce convergent estimates of surface fluxes with increasing grid sizes. Our finding demonstrates that an atmospheric model can underestimate surface fluxes in regional scale not necessarily due to intrinsic model inaccuracy (e.g., inaccurate parameterization) but due to scale-dependent nonlinear effect of spatial variability in atmospheric conditions.


Introduction
Surface fluxes from the earth surface are playing a key role in climate changes by modifying the circulation of the overlying atmosphere (e.g., Avissar et al., 1989;Baldocchi et al., 2001).Regional networks of direct surface flux measurements have helped significantly to improve our understanding of soil-vegetation-atmosphere (SVAT) interactions (e.g., Baldocchi et al., 2001).In particular, they have provided a pivotal information in develping and evaluating SVAT models at plot scale (≤1 km 2 ).So far, continuous estimation of surface fluxes is not robust beyond plot scale (1-10 km 2 ) and our understanding of the SVAT interactions remains a major Correspondence to: Jinkyu Hong (jkhong@yonsei.kr)challenge in the regional scale (Dolman et al., 2006).Atmospheric mesoscale models have been a powerful tool and have a potential to monitor energy and mass exchanges in regional scale although modeling results are yet to be validated (e.g., Henderson-Sellers et al., 1993;Oncley and Dudhia, 1995;Chen and Dudhia, 2001a;Kerschgens and Heinemann, 2005; van der Molen and Dolman, 2007).
When we interpret the simulated SVAT processes in regional scale, it is requisite to scale down the model results because of the lack of the observation data to cover regional scale.Also, we should note that most of atmospheric numerical models implicitly assume spatial homogeneity of surface properties as well as atmospheric conditions such that they calculate gridbox averaged surface fluxes using mean variables of surface properties and meteorological fields.Indeed, spatial heterogeneity in surface properties (e.g., vegetation types) (Sellers et al., 1995;Middelburg et al., 1999;Baldocchi et al., 2005;Li et al., 2008) and in atmospheric conditions (e.g., radiative fluxes) has been a critical issue in validating numerical models and remote sensing (Bertoldi et al., 2008).
With respect to this impact of spatial heterogeneity on the model evaluation, two aspects of scaling issue are worth noting.One is related to a small size of tower footprint compared to a typical grid size of a mesoscale model.Such a difference can induce inadvertent model bias despite proper subgrid parameterizations in the model (Oncley and Dudhia, 1995;Kim et al., 2006).Previous studies, however, validated the models with grid sizes from 10 to 200 km against the insitu tower data of ∼1 km footprint without clear assessment of surface heterogeneity around a measurement location and in the model domain (e.g., Chen and Dudhia, 2001a;Brotzge, 2004;Pyles et al., 2003).
The other issue is the nonlinear effect due to spatial variability in surface properties and/or atmospheric conditions expressed by Jensen's inequality in mathematics (Fig. 1) (Sellers et al., 1992(Sellers et al., , 1995)).Following the pioneering work of Sellers et al. (1992), many models now consider at Published by Copernicus Publications on behalf of the European Geosciences Union and the American Geophysical Union.rface fluxes, that of atmospheric conditions received little ention.Similar to the effect of surface heterogeneity on the ulated surface fluxes, spatial variability in meteorological lds for driving SVAT modules in a fully coupled mesoscale odel can also result in potential bias of the modeled surface xes with different model grid size.To properly assess surface fluxes at regional scales using ully coupled atmospheric model, we attempted to quantify nlinear effect due to the spatial variability in atmospheric nditions.Without consideration of this scaling issue, exination of the model performance on surface fluxes could d us to a hasty conclusion that a model output agrees well th field observation or that the model has systematic bias subgrid parameterization.As far as we know, 'scale-invariance' of surface fluxes in lation to spatial variability in atmospheric conditions has t been tested in fully coupled mesoscale models.Here, the ale-invariance of surface fluxes implies that surface fluxes gregated from a fine grid should be equal to those from a arse grid (Fig. 2).Therefore, such a test offers the inforation to constrain the simulated surface fluxes.Based on e assessment of the scale-dependency, we have evaluated e SVAT processes in the model at a tower footprint scale 2) (see Sellers et al. (1992Sellers et al. ( , 1997) ) for tion).If F is surface flux as a functio surface conductance, soil moisture), Thus, complete values of F over an lated from where A is the area of domain and t average over A. However, because of 115 an entire domain is calculated: where ⟨F ⟩ M is an estimate of ⟨F ⟩.
Ideally, surface fluxes aggregate should be consistent with values from 120 and vegetation have, however, non least surface heterogeneity by combining different vegetation types within a grid (e.g., Bonan et al., 2002).However, compared to popular interests in the impact of spatial heterogeneity of surface properties (e.g., vegetation type) on the simulated surface fluxes, that of atmospheric conditions received little attention.Similar to the effect of surface heterogeneity on the simulated surface fluxes, spatial variability in meteorological fields for driving SVAT modules in a fully coupled mesoscale model can also result in potential bias of the modeled surface fluxes with different model grid size.
To properly assess surface fluxes at regional scales using a fully coupled atmospheric model, we attempted to quantify nonlinear effect due to the spatial variability in atmospheric conditions.Without consideration of this scaling issue, examination of the model performance on surface fluxes could lead us to a hasty conclusion that a model output agrees well with field observation or that the model has systematic bias in subgrid parameterization.
As far as we know, "scale-invariance" of surface fluxes in relation to spatial variability in atmospheric conditions has not been tested in fully coupled mesoscale models.Here, the scale-invariance of surface fluxes implies that surface fluxes aggregated from a fine grid should be equal to those from a coarse grid (Fig. 2).Therefore, such a test offers the information to constrain the simulated surface fluxes.Based on the assessment of the scale-dependency, we have evaluated the SVAT processes in the model at a tower footprint scale using in-situ tower data.

Theoretical background: scale-invariance of surface fluxes
The concept of scale-invariance of surface fluxes was first proposed by Sellers et al. (1992).The test of scale-invariance can give information to constrain surface fluxes in the model.The issue of scale-invariance with respect to model performance is the extent to which the calculated surface energy balance (SEB) at coarse domain is influenced by the use of coarse-domain averaged surface boundary conditions (Fig. 2) (see Sellers et al. (1992Sellers et al. ( , 1997) ) for more detailed information).If F is surface flux as a function of x i (e.g.topography, surface conductance, soil moisture), we can write Thus, complete values of F over an entire domain are calculated from where A is the area of domain and the brackets denote area average over A. However, because of practical reason, F over an entire domain is calculated: where F M is an estimate of F .Ideally, surface fluxes aggregated from fine grid size should be consistent with values from coarser grid size.Soil and vegetation have, however, nonlinear responses to the atmosphere and consequently surface fluxes are not conserved (i.e., F = F M ).Consequently, testing the scaleinvariance of surface fluxes to the degree of spatial variability in atmospheric conditions can provide valuable information on the model evaluation.

Field observations
We used the data collected at the ANNI flux station near Naqu, Tibet (31.37 • N; 91.90 • E, 4580 m above m.s.l.) from 26 to 31 August 2002.The site was a flat and homogeneous shortgrass prairie (Fig. 3).The soil was predominantly sandy silt loam and soil surface is sparsely covered with short grass (canopy height of <0.05 m and leaf area index of <0.5) due to grazing.
We used the eddy-covariance system mounted on a 20 m tower to measure the turbulent fluxes of sensible heat (H) and latent heat (LE).The system consisted of three-dimensional sonic anemometer (CSAT3, Campbell Scientific Inc., US) and infrared gas analyzer (LI7500, LiCor, US).Sampling rate was 10 Hz and averaging time was 30 min.Other supporting measurements included surface radiative fluxes, soil moisture, soil temperature, and precipitation.During the study period, energy balance was nearly closed and quality assurance has been done based on Lee et al. (2004).More information on field observations at the ANNI flux station can be found in Hong et al. (2004), Hong and Kim (2008) and http://ceop.net.

Meteorological conditions
During the simulation period, the site was under the high pressure system centered on the eastern Tibet.The weather charts show that the position and intensity of this high pressure system had fluctuated daily and this system over the Plateau scaled down from 28 August 2002.Subsequently, the atmospheric surface pressure tended to decrease at the flux station.The atmospheric conditions at the flux station were similar to the weather chart reports and NCEP-NCAR reanalysis data (hereafter NNRD) except rainfall.Based on the weather charts, it was overcast on 28 and 29 August but there was no precipitation on the Plateau.However, 0.6 mm rainfall was recorded for an hour at the flux station in late evening of 28 August.After this rainfall, the mixing ratio increased and then slowly decreased with time and the Bowen ratio kept increasing from about 0.3 to 0.5 with decreasing soil moisture.

Numerical experiment
In this study, for testing scale-invariance in the atmospheric mesoscale model, we carried out the numerical simulation using Penn State/NCAR mesoscale model, MM5 from 26 to 31 August 2002.Using a nesting option, the model was run on a 60 km horizontal grid (domain 1) with subdomains of 20 and 1 km horizontal grid increments (domain 2 and domain 3, respectively) (Fig. 4).To assess the scale dependence of surface fluxes (Eq.3), the outputs from the domain 3 were aggregated to 20 km grid scale.These aggregated outputs were compared with the outputs from the domain 2 simulation.
The number of vertical layers was 35 and the model top was located at 50 hPa.The initial and boundary conditions were based on NNRD.Topography data were 30 s tiled topography data from USGS (US Geological Survey) and monthly vegetation fraction data were from AVHRR (Advanced Very High Resolution Radiometer) (Gutman and Ignatov, 1998).For land use and soil classification, 25category USGS vegetation data and STATSGO (State Soil Geographic) 17-category soil data were utilized.The Kain-Fritsch cumulus parameterization scheme (Kain and Fritsch, 1993) and Reisner2 microphysics scheme were used (Reisner et al., 1998) in this study.But cumulus parameterization was turned off in domain 3 because 1 km horizontal grid spacing in domain 3 is assumed to be small enough to resolve the cumulus convection.We used MRF PBL scheme and RRTM radiation scheme (Hong and Pan, 1996 , 1997).For land surface processes, Noah land surface model (LSM) (Koren et al., 1999) was used in this study.More information can be found on the MM5 website (http://www.mmm.ucar.edu/mm5/).The measurement height was 20 m for tower observation but model outputs were from 2 m for temperature and humidity and 10 m for wind speed.For proper comparisons, the observation data were recalculated using Monin-Obukhov theory (Kaimal and Finnigan, 1994).The boundary layer processes used in this model are based on Monin-Obukhov (MO) similarity and assume that source/sink distribution is same among different scalars (i.e., temperature and water).Recently, McNaughton and Brunet (2002) suggested that surface fluxes can deviate from the prediction of MO similarity.This issue was observationally investigated by Hong et al. (2004) and Choi et al. (2004) but such a deviation does not impact on our conclusion in this study.

Footprint analysis
Turbulent flux and scalar concentration measurements are influenced by the underlying surface below the sensors and represent a weighted average of upwind conditions.A measured quantity, η, is given to the integrated sum of source strength distribution (Q n ) at location r and footprint (or source weight distribution), f (Schmid, 1997): The footprint (f ) determines the relative weight of individual point sources.To obtain the spatial context affecting the measured quantities, the source area of level P is given as: where the source area ( P ) is the surface area bounded by a footprint isopleth f (x, y, z m ), such that P is the fraction of the total integrated footprint function (φ tot ) and φ p is the integral of the footprint over P (Schmid, 1997).The underlying surface area of influence of sensors continuously changes with measurement height, atmospheric stability, wind direction, and turbulence structure.
In this study, we calculated a tower footprint using fluxsource area model (FSAM) by Schmid (1997).The outputs of FSAM are the characteristic dimensions of the footprint such as maximum source location (x m ), near end (a), far end (e) and maximum lateral half-width (d) of the source area (Fig. 5).During the simulation period, slightly unstable conditions occurred more frequently and data were evenly distributed for the wind directions except northeasterly.During the daytime, average e and d were about 400 m and 40 m.On the other hand, in the nighttime, average e and d increased up to about 7 km and 4 km, respectively.

Test of scale-invariance of surface fluxes
Figure 6 presents time series of sensible (H ) and latent heat fluxes (LE) from the tower observation, reanalysis data and the model.H and LE from domain 3 (1 km grid size) reproduced the observed patterns and their magnitudes compared relatively well to those from domain 1 (20 km grid size) (Table 1).However, the overestimation of cloud and precipitation amount in the model caused the poor model performance, and yet the inappropriately simulated pattern in domain 1 produced a relatively good performance in rainy days.Here, we note that the surface fluxes at domain 1 and 3 were comparable in the beginning of the simulation and then diverged with time.This emphasizes that the model outputs can have inconsistent patterns with different grid sizes compared to the observation.
This discrepancy of the outputs from different grids was proportional to the domain averaged standard deviation of surface fluxes in the model (Table 2).The domain averaged standard deviation increased as the numerical integration progressed, indicating the enhanced spatial variability with time.Two-dimensional Fourier transform also corroborates that the spatial variabilities of surface fluxes increased in the model as the numerical integration progressed (Fig. 8).We speculate that this is the result of homogeneous initial and boundary conditions disaggregated down from the coarse domain.In the early period of the simulation, initial and boundary conditions for domain 3 were obtained from domain 2 by linear interpolation and thus had relatively homogenous field.As the numerical integration progressed, the model eventually made its own heterogeneity of meteorological fields and surface fluxes.It also suggests that a spin-up process is Nonlin.Processes Geophys., 15,[965][966][967][968][969][970][971][972][973][974][975]2008 www.nonlin-processes-geophys.net/15/965/2008/ required at least for one day so that the model adapts itself from the smooth boundary conditions.
Our analysis shows that this increased spatial variability of H and LE is directly related to spatial distribution of clouds and precipitation in the model.The amount of cloud and precipitation directly affects the downward radiation at the surface, which in turn limits the available energy for other surface fluxes such as H and LE (Fig. 7).Until now it was not certain that such an increased heterogeneity reflects the real environment because we only had the data to cover about 1 km domain.However, the domain 3 outputs well captured the observed patterns and magnitudes of surface fluxes and this increased heterogeneity eventually resulted in the failure of the scale-invariance of surface fluxes in the larger domain (Table 3).
A plausible explanation of this failure can be explained the surface radiatve forcing due to the spatial variability in clouds and rainfall in the model.Surface fluxes respond to radiative forcing as the cosine function (Sellers et al., 1997).In case of spatial variability of radiative forcing caused by moderate terrain, Sellers et al. (1997) showed that this nonlinear effect can be negligible.However, in our study, surface radiation spatially spread over wide ranges because of the enhanced spatial variability in clouds and precipitation, and consequently the nonlinearity was not negligible.Under ).Here M i and O i are the simulated and observed values respectively.d is 1 for a perfect model (Willmott, 1982) this condition, we can no longer approximate the response of SEB to such varying radiation with a linear function and thus the scale-invariance is not satisfied.Also, unlike Sellers et al. (1997), spatial variability of soil water contents did not decrease as the study area dried out because of the local scale precipitation in the domain and the patchiness of LE (Fig. 7).Under these conditions, the model does not provide convergent estimate of surface fluxes in coarse grid simulation.

Model validation using a tower observation
Overall, the domain 3 simulation credibly reproduced the magnitudes and patterns of surface fluxes on clear days (Table 1).The model, however, did not properly simulate precipitation and clouds even in domain 3 (Table 1; Fig. 9).On 28 August, there were afternoon rainfalls in both the observation and the model and the model captured well the immediate reduction of the downward shortwave radiation due to this rainfall.However, the total amount of precipitation on 28 August was different from the observed precipitation.Furthermore, the duration of the rainfall in the model was longer (∼5 h) than that of the observation (∼1 h) and larger amount of clouds continued for longer period in the model.On 27, 29, and 30 August, there were small amount of rainfalls in the model (<2.5 mm/30 min) but these rainfalls were not observed at the flux station.Such bias of precipitation in the model substantially reduced the downward radiative fluxes and other surface fluxes through the reduction of available energy (  The rainfall in domain 3 (1 km resolution) around the tower was the grid resolvable rainfall, but not the grid resolvable rainfall in domain 2 (20 km resolution).It indicates that the convective system was of the order of a few kilometers in the model.On the other hand both convective and non-convective rainfall existed together in domain 1 (60 km resolution), suggesting that precipitation had different origin and features with different grid sizes.
While previous studies have reported the overestimation of downward shortwave radiation in mesoscale models, which was attributed to the lack of aerosol parameterization in the radiation schemes (e.g., Betts et al., 1997;Chen and Dudhia, 2001b), the simulated downward shortwave radiation was consistently larger except around noon on clear days.Our result indicates that the optical depth in the model was shallower at high solar zenith angle but thicker at low solar zenith angle over the high elevation area like the study site (Fig. 9).Indeed, Reiter et al. (1987) reported that the atmosphere over the study site was highly transmissive and there was an additional scattering source from cumulus clouds, thus overestimating the diffusive radiation at high solar zenith angle.
Similarly, the model also overestimated the downward longwave radiation, indicating that the modeled atmosphere was warmer than the observation.In particular, such an overestimation was more pronounced when there was rainfall.Chen and Dudhia (2001b) also reported a similar overestimation of downward longwave radiation in the MM5 mesoscale model.In the radiative transfer scheme used in this study, the downward longwave radiation depends on the amount of water vapor in the atmospheric column through the modification of emissivity.Therefore, we speculate that such bias of downward longwave radiation was the result of the overestimated clouds and water vapor in the model.
The simulated wind speeds were smaller than the observation in general (not shown here).On the contrary the model overestimated friction velocity (u * ) except for the underestimation in the nighttime conditions due to the minimum u * assigned in Noah LSM.This manifests that roughness length in the model was relatively large, thereby reducing wind speed excessively.Also, this suggests that the minimum friction velocity in the Noah LSM should be modified.Interestingly, the relative deviation of the simulated u * from the observed u * showed the scale-dependence on the tower footprint (Fig. 10).Tower footprints were >1 km in nighttime and probably include more heterogeneous surface than well mixed daytime turbulence conditions.Accordingly, such dependence might be related to the scaling issue due to surface heterogeneity or the model's systematic bias between day and night.To accurately address this scaling issue, the potential bias of flux observations due to the footprint mismatch should be quantified with satellite image analysis (Kim et al., 2006).
Figure 11 shows diurnal variations of surface temperature and soil water content in the model with the observation.Similar to the simulated downward shortwave radiation, the model produced surface skin temperature (T s ) on clear days comparable to to the observation but T s decreased rapidly near sunset due to the inappropriate rainfall simulation.In case of soil moisture, the model systematically underestimated soil water content.The observed soil water content decreased monotonically with time during the study period.The soil water content in the model also decreased with time but slightly increased due to the precipitation on 28 August.The initial value of soil water content in the model deviated from the observed values but this discrepancy was reduced mainly due to the excessive rainfall in the model.It is also worth noting that the model reproduced LE relatively well despite the systematically smaller soil water contents than the observation.This is likely through strong atmospheric evaporative demands in the model.Indeed the simulated relative humidity was smaller than the observed value, and the reduction of soil water content for LE was partially compensated with the excessive rainfall in the model.

Summary and Conclusions
We investigated the scale-dependency of surface fluxes due to spatial variability of atmospheric conditions in the mesoscale model.To properly assess the scale-dependency with different grid sizes, we performed a numerical simulation from a tower footprint to regional scale using the atmospheric mesoscale model.
Overall, the mesoscale model well captured the diurnal variation of surface fluxes on clear days in the spatial scales of the tower footprint.However, large scale atmospheric features, which produced the conditions favorable to convective systems in the coarse domain, induced the overestimation of the amount of clouds and precipitation even in the fine resolution.Eventually the growth of those convective systems led to the increase of spatial heterogeneity in surface radiative forcing within the model domain.As the degree of spatial heterogeneity in the meteorological fields grew up in the model, surface fluxes from different grid sizes showed the scale-dependency.Accordingly, the model did not provide convergent estimates of surface fluxes to the aggregation and disaggregation processes in different grid sizes.Sellers et al. (1997) showed the insignificant impacts of surface heterogeneity on surface fluxes in the off-line biosphere model by assuming that atmospheric conditions were relatively homogeneous.In our study, however, substantial spatial variability in atmospheric conditions was generated in the fully coupled mesoscale model, which ultimately led to the failure of scale-invariance.Furthermore, our study demonstrates that the surface heterogeneity in soil moisture and temperature was imposed by larger scale constraints such as the convective rain and radiative forcing in the model.
Nonlin.Processes Geophys., 15,[965][966][967][968][969][970][971][972][973][974][975]2008 www.nonlin-processes-geophys.net/15/965/2008/ Our findings have an important implication when we evaluate surface fluxes in regional and continental scales using the atmospheric model: 1) One should carefully apply the modeling outputs from a coarse domain for inferring regional scale surface fluxes because the modeled surface fluxes can be underestimated due to the failure of scale-invariance as the spatial heterogeneity grows up in the model; and 2) any difference between model outputs and the observation data can result from the failure of scale-invariance.That is, the model bias will not necessarily result from the intrinsic model inaccuracy when spatial variability in atmospheric conditions prevails.Only when the scale-invariance requirement is satisfied, one can relate the model biases to the structural deficiencies of model.In case of the scale-dependent surface fluxes, we may need to add the scale-dependent areaaveraged forms of the parameterization for a proper assessment of regional surface fluxes (Sellers et al., 1997).Otherwise, we need to aggregate surface fluxes from the fine resolution simulation of the tower footprint.
Jensen's inequality.F is a function of x and <> is an eraging operator.

Fig. 2 .
Fig. 2. The concept of (dis)aggregation a face fluxes in the model adapted from Se

Fig. 1 .
Fig. 1.Jensen's inequality.F is a function of x and <> is an averaging operator.

Fig. 2 .
Fig. 2. The concept of (dis)aggregation and scale-invariance of surface fluxes in the model adapted from Sellers et al. (1992).

Fig. 3 .
Fig. 3. Eddy-covariance system was installed at 20 m tower in the ANNI flux station.The site was a flat and homogeneous shortgrass prairie.

Fig. 5 .
Fig. 5. Characteristic dimensions of the source area of FSAM: maximum source location (x m ), near end (a), far end (e), and maximum lateral half-width (d) of the source area.

Fig. 8 .
Fig. 8. Base-10 logarithm of 2-dimensional power spectrum density of LE and H at 05:00 UTC on 26 and 30 August.k x and k y is a wave number of longitudinal and latitudinal direction, respectively.

Fig. 9 .
Fig. 9. Comparison of the simulated (domain 3) and observed downward longwave and shortwave radiation and precipitation.

Fig. 10 .
Fig. 10.Ratio of u * differences between the model and observation with far end distance, e.

Table 1 .
Index of agreement d of the simulated surface fluxes (≡1− (

Table 2 .
. Standard deviation of LE (σ LE ) and H (σ H ) from the domain 3 averaged mean values at 05:00 UTC.

Table 3 .
Comparison of LE, H, downward shortwave radiation (R sdn ) and soil water content in domain 2 (20 km grid size) and in domain 1 after aggregating to 20 km spatial scale around a flux tower at 05:00 UTC.Value in a parenthesis is a standard error from the spatial average.D2, D3, and Obs are the data from domain2, domain3, and observation, respectively.