An Energy Balance Model for Low--Level Clouds Based on a Simulation Resolving Mesoscale Motions

This study proposes a new energy balance model to determine the cloud fraction of low-level clouds. It is assumed that the horizontal cloud field consists of several individual cloud cells with a similar structure. Using a high–resolution simulation dataset with a wide numerical domain, we conducted an energy budget analysis. Here we show that the energy injected into the domain by surface flux is approximately balanced with that loss due to radiation and advection due to large–scale motion. The analysis of cloud cells within the simulated cloud field showed that the cloud field consists of a number of cloud cells with similar structures. We developed a simple model for the cloud fraction from the energy conservation equation. The cloud fraction determined using the model developed in this study was able to quantitatively captured the simulated cloud fraction.


Introduction
Low-level clouds, such as stratocumulus or shallow cumulus, play an important role in determining the radiation budget of the globe because such clouds cover a wide area and reflect shortwave radiation (Klein and Hartmann 1993;Klein et al. 1995;Wood 2012;Zelinka et al. 2016). These clouds exist at low altitudes (z < 2 km) with a depth of several hundred meters. Satellite observations have revealed that the horizontal fields of these low-level clouds have cellular structures and that there are several kinds of horizontal cloud patterns (Krueger and Fritz 1961;Stevens et al. 2005b;Wood and Hartmann 2006;Comstock et al. 2007). Two major patterns are open cell and closed cell structures (Koren and Feingold 2013). The former pattern accompanies a cumulus-like flow and cloud structures, i.e., narrow updraft and wide downdraft regions. Horizontally, the low-level cloud field covers hundreds of square kilometers and consists of a number of cloudy cells in which each cell has a circulation with a horizontal scale of several kilometers. The large areas of the planet covered by these low-level clouds have a significant influence on the global radiation budget.
Various physical processes that cover a wide range of spatiotemporal scales are important for the development and maintenance of low-level clouds, and these range from interaction between aerosol and cloud droplets to large-scale motion (e.g., the subsidence branch of the Hadley circulation). (Bretherton and Wyant 1997;Wyant et al. 1997;Bretherton et al. 1999) For instance, when the number of cloud condensation nuclei (CCN) is small, large water droplets form in a saturated environment, and these are more likely to fall as rain. This precipitation (drizzle) greatly changes the cloud structure and causes a structural transition from closed cell to open cell (Xue et al. 2008;Caldwell and Bretherton 2008;Stevens and Feingold 2009;Wang and Feingold 2009a, b;Feingold et al. 2010Feingold et al. , 2015Berner et al. 2011;Yamaguchi and Feingold 2015). Another example representing the importance of various scales for low-level clouds is buoyancy reversal near a cloud top, in which entrainment of dry air aloft into the cloud layer promotes evaporation and cooling, and downward motion is accelerated, resulting in a reducing stratocumulus cloud deck (Yamaguchi and Randall 2008;Mellado 2010;Mellado et al. 2010Mellado et al. , 2014Noda et al. 2013van der Dussen et al. 2014).
Because of the importance of small-scale microphysical processes, it is difficult to explicitly resolve the low-level clouds in global atmospheric models. Even a global simulation with a grid spacing of 870 m (Miyamoto et al. 2013;Kajikawa et al. 2016;Yashiro et al. 2016) has difficulty to accurately simulate the low-level clouds and hence cloud fraction. A number of studies have investigated low-level clouds using the grids with a resolution of O(10) m (Moeng et al. 1996;Stevens et al. 2005a;Ackerman et al. 2009;Yamagu-chi and Feingold 2012;Sato et al. 2015a, b). These studies successfully simulated the low-level clouds, i.e., the resolved scale needs to be in the order of 10 m (unresolved scale is less than this scale) to realistically simulate the low-level clouds. In other words, one of the dominant mechanisms governing these clouds is of this order. As the global radiation budget is affected by low-level clouds, global models, especially for long-term simulations such as climate simulations, need to properly incorporate their effects, and hence a deeper understanding of these low-level clouds would be beneficial.
Previous studies have studied energy and water budget in the atmospheric boundary layer by assuming an equilibrium state (Lilly 1968;Schubert 1976;Albrecht et al. 1979;Betts 1976;Betts and Ridgway 1989;Neggers et al. 2006). Caldwell et al. (2005) conducted budget analyses for mass, heat and liquid water static energy in the boundary layer by observation and reanalysis data to examine the entrainment at the top of boundary layer. Kalmus et al. (2014) also conducted budget analyses based on a set of satellite, GPS, and ship-based data. They found that in climatological mean, the transition from stratocumulus to cumulus state is associated with an increase in surface latent heat flux, boundary layer height, rain, and horizontal advection of dry air and a decrease in entrainment of dry air. Chung et al. (2012) showed from an energy budget analysis based on a series of large-eddy simulations (LESs) that the temperature tendency due to radiative cooling is balanced with the tendency due to the subsidence warming during the transition of cloud regime. They derived an equation for cloud fraction based on the balance of two tendency terms. However, the numerical domain of their simulation covered horizontally 3.2 ´ 3.2 km 2 , which is not large enough an area to resolve mesoscale motions.
A notable feature of low-level clouds is that the cloud field consists of a number of cloudy cells, each of which has similar flow and cloud structures. The open-and closed cell flow fields are similar to the Rayleigh-Bénard convection (Krishnamurti 1975;Laufersweiler and Shirer 1989;Weidauer et al. 2010Weidauer et al. , 2011Miyamoto et al. 2020). The convection transports heat vertically via a number of cloud cells that have a flow structure similar to each other. In fact, many of the previous studies on low clouds introduced above implicitly state that the cloud field consists of similar cloud cells. This unique feature of low-level clouds is one of the key assumptions on which the present study is based.
In this study, we develop an energy balance model for low-level clouds, based on the key assumption that the cloud field consists of a number of cloud cells with the same structure, and we use the model to conduct an energy budget analysis of a simulated cloud field. We use the high-resolution simulation with a wide numerical domain conducted by Sato et al. (2015b) for our analysis. The simulation setting and methodology used to analyze the cloud cells are described in Section 2. An overview of the simulation of Sato et al. (2015b) and the results of the detected cloud cells are presented in Section 3. An energy budget analysis is performed and an energy balance model is presented in Section 4. The results are discussed in Section 5 and we present our conclusions in Section 6. Sato et al. (2015b) We analyzed the results of the idealized numerical simulation in Sato et al. (2015b), which covers (768, 28, 2) km in the (x*, y, z) directions with grid intervals of (50, 50, 5) m. A fully compressible numerical model, Scalable Computing for Advanced Library and Environment (SCALE) (Nishizawa et al. 2015;Sato et al. 2015a), was used for the integration. The prognostic quantities generated by SCALE were the density of the total mass, three-dimensional momentum, potential temperature weighted by density, and microphysical quantities (mixing ratio of water vapor, and mass and number density of cloud water, rain water, ice, snow, and graupel). The time differential was discretized using the three-step Runge-Kutta scheme. The advection and pressure gradient terms were discretized using the fourth-and second-order accuracies, respectively. The discretized equations were solved explicitly for both the horizontal and vertical directions.

Simulation design of
The number of grid squares was 6144 ´ 564 ´ 275. The grid spacing was vertically stretched above 1.2 km, and also horizontally in 0 km < x* < 247 km, and 545 km < x* < 768 km. From x* = 247 km to 545 km, the grids were evenly allocated every 50 m, which is the analysis domain and hereafter we focus on this region. x* = 245 km is defined as the upwind boundary (x = 0). The boundary conditions for the x* and y directions were open and periodic, respectively. The effects of sub-grid-scale turbulence were solved using a Smagorinsky scheme generalized for anisotropic grids (Smagorinsky 1963;Lilly 1962;Scotti et al. 1993). Cloud physical processes were calculated using double moment-bulk cloud microphysics (Seifert and Beheng 2006;Seiki and Nakajima 2014). The nucle-ation was represented by assuming a temporally constant number of CCN. Longwave radiation was only considered in the simulation and vertical radiation fluxes were solved following Stevens et al. (2005a). The inversion height used in the radiation scheme was defined as the level at which the total water mixing ratio q t became less than 8.0 g kg −1 . In the horizontally stretched regions and the topmost 500-m depth of the domain, Rayleigh damping was applied to all of the prognostic variables to prevent artificial reflection of gravity waves. The timescales used for the damping were 300 s and 10 s for the horizontal and vertical directions, respectively.
The initial vertical profiles for the temperature T and water vapor mixing ratio q v were constructed from an observation campaign, Second DYnamics and Chemistry Of Marine Stratocumulus (DYCOMS-II) RF02 (Ackerman et al. 2009). The effect by largescale subsidence (w LS = Dz) was given to all prognostic variables (Ackerman et al. 2009), where D was 1.33 ´ 10 −6 s −1 , the same as that used in Berner et al. (2011). The integration period was 16 hours. The numerical simulation was initialized from horizontally uniform fields for all quantities except the surface heat fluxes and number of CCN, which were changed in the streamwise (x*) direction. Surface fluxes for sensible and latent heat were 15 W m −2 and 93 W m −2 at x* = 247 km, and the fluxes increased at a rate of 0.03062 W m −2 km −1 and 0.1365 W m −2 km −1 in the x* direction, respectively. Thus, the sensible and latent heat fluxes were 24.12476 W m −2 and 133.677 W m −2 at the downwind edge. The equations for surface fluxes for sensible and latent heat were given by and hence the surface enthalpy flux was given by CCN decreased according to 250 exp(-7.0433x* ´ 10 −6 ) cm −3 and CCN at the downwind edges was 5.38094 cm −3 . Both the surface fluxes and CCN were fixed during the simulation. A uniform velocity of 5 m s −1 was present in the x* direction in the initial field. The parameterization of Stevens et al. (2005a) was used for the net radiative longwave flux.
Here, Q 1 (z) = 85 ò z ztop ρq l dz¢ and Q 2 (z) = 85 ò 0 z ρq l dz¢ are the integrated liquid water mixing ratio in the vertical direction, F R0 = 70 W m −2 , F R1 = 22 W m −2 , ρ z i is the air density at the cloud top level z i , D s = 3.75 ´ 10 −6 s −1 , and a z = 1 m −4/3 . More specific information on the experimental setting can be found in (Sato et al. 2015b).

Detection algorithm of for cloud cells
We developed a method that is able to detect cloud cells in a simulated cloud field based on an approach designed for deep convection (Miyamoto et al. 2013. First, vertical velocity was vertically averaged in the boundary layer, and horizontally smoothed 100 times by applying a 1-2-1 filter. Second, the grid point of the cell center was defined as a grid point having a local peak of vertical velocity relative to the standard deviation of vertical velocity obtained in a subdomain, σ w . Specifically, grid points, at which the absolute vertical velocity exceeded the standard deviation; i.e., |w | > σ w , were detected as the cell center. We used the standard deviation to search for strongly (anomalously) deviated peaks in the subdomains. Since the methodology applies the absolute velocity, the method can capture positive and negative peaks in vertical velocity. In the simulation, positive peaks were detected. Once a center grid was detected, the coordinates were transformed into cylindrical coordinates around the detected cell center. Sato et al. (2015) Figure 1 shows the horizontal distribution of the liquid water path (LWP) at seven selected simulation times. The cloud field shows cellular structures at all times. The spatial scale, or distance between cellular structures, increases with time, especially on the downwind side (right hand side of panels). At the same time the spatial scale of the cloud cell itself becomes also larger. The horizontal cloud field on the upwind side looks like the cloud structure of a closed The vertical velocity is positively large at the edge of the cellular structure, especially in R03-R07. Figure 3 shows a time series of the cloud fraction, which is defined as the fractional area where LWP is greater than 80 g m −2 in each subdomain. The cloud fraction is close to 1.0 in all subdomains at the beginning of the simulation, which indicates that the entire domain is covered by clouds immediately after the simulation begins. The temporal changes in cloud fraction are similar in all subdomains, whereas the magnitudes differ; i.e., the magnitude decreases after integration begins and then remains at a constant value after t = 12 h. On the upwind side (R01 or R02), the magnitude remains large throughout the simulation. In contrast, the cloud fraction rapidly decreases to 0.7 on the downwind side (R06-R09). Figure 4 shows vertical profiles of the temperature, water vapor mixing ratio, cloud water mixing ratio q c , rain water mixing ratio q r , and vertical fluxes of the four quantities in R01, R03, R05, and R07. The quantities are temporally averaged from t = 12 h to 16 h and horizontally averaged, whereas the fluxes are integrated in each subdomain. T decreases with height up to z = 0.85 km, then rapidly increases to 292 K, and is constant above. q v has a peak at the surface, decreases with height, rapidly decreases at z = 0.85 km, and is nearly constant above. q v is large, especially at lower levels, which most likely results from surface fluxes. The larger T and q v in the downwind region are caused by the increased surface heat flux. q c has a peak below the inversion height and it is large in the upwind region. The vertical profiles of q r are similar in the four subdomains, showing a peak below that of q c and a non-zero value at the surface.

Overview of
Vertical profiles of the area-integrated vertical flux of the four quantities in the four subdomains are shown in Figs  q c is largest at the peak altitude of q c , and is large in the upwind region. The vertical profi le of the fl ux of q r has a negative peak at z = 0.28 km and a positive peak at z = 0.78 km. We note that the fl ux of q r does not include the precipitation fl ux in which the water droplets fall by gravity.

Structure of cloud cells
In order to examine the cell structure in the subdomains, we developed a methodology for cell detection. The cloud cell detection method was applied to the vertical velocity fi eld (cf. Fig. 2) from t = 12 h to 16 h. Figure 5 shows the number of detected cells for each region averaged over the analysis period in each subdomain. The number is larger on the upwind side, with the largest being 12 in R02, whereas it is 7 -8 on the downwind side. The tendency of the number of cells is consistent with a visual inspection of Figs. 1 and 2.
To evaluate the proposed method, we tested the sensitivity of the number of detected cells by changing the number of smoothing processes and the threshold for the vertical velocity. In particular, by introducing a factor, B, to the defi nition |w | > Bσ w , two values, B = 0.5 and 2.0, were tested. The top panels in Fig.  6 show the number of detected cells when B = 0.5 (left) and B = 2.0 (right). Overall, the number of cells was not largely changed by the thresholds, whereas the sensitivities of the number of detected cells to the two parameters were reasonable. Compared with the control value B = 1.0 (Fig. 5), the number slightly increased and decreased when the factor B decreased . Vertical profi les of (upper) temperature (K), water vapor mixing ratio (g kg −1 ), cloud water mixing ratio (g kg −1 ), rain water mixing ratio (g kg −1 ), liquid-water potential temperature (K), and (lower) vertical fl uxes of temperature (K m s −1 ) and water mixing ratios (g kg −1 m s −1 ). The quantities in the upper panels are horizontally averaged and the fl uxes in the lower panels are horizontally integrated in each subdomain. They are averaged from t = 12 h to 16 h. The dashed lines in the panels, (a), (b), and (e), indicate the initial values. and increased, respectively. Next, the number of smoothing N was also tested by applying N = 50 and 200, of which results of detection are shown in the bottom panels in Fig. 6. The number of cells was not largely sensitive to the number of smoothing; the cell number was slightly reduced when the number of smoothing was 200, while the number did not change when the number of smoothing was reduced to 50 times. Figure 7 shows a radius-height cross section of the radial and vertical velocities, which are composites of all of the cloud cells detected. The radial velocity has a positive peak immediately below the inversion height, and a negative peak above the surface, both of which are located at a radius of several hundred meters. The vertical velocity has a positive peak at the cell center.
The outside is the downward region, while the magnitude of the downward velocity is much smaller than that of upward velocity. The flow fields in the vertical cross section are qualitatively the same in the subdomains. This indicates that circulation is present in the cell with the inward flow at the bottom, upward motion at the convection center, and outward motion at the top of the boundary layer. The flow field is confined to a small area on the upwind side, whereas the magnitude of the vertical velocity is almost the same among the subdomains.
Radius-height sections of q c and q r are shown in Fig. 8. The detected cells in all subdomains have a cloud layer approximately from z = 550 m to 850 m, a peak of w at the core, and q r immediately outside the core. The region of high q r extends from the upper layer to the surface in the downwind region, indicating that the precipitation reaches the ground. The precipitation is large on the downwind side. The composites of velocity and water substance have also been produced by Zhou and Bretherton (2019). The present result is qualitatively the same as the previous study, whereas the methodology to construct the composite and the simulation data are different. Figure 9 shows radius-height cross sections of the number densities of cloud water N c and rain water N r . A large number density is approximately collocated with the high mixing ratio, whereas the peak of N c  is located lower than that of q c . They have a peak at the cell center that decays radially. In particular, N c is large in the upwind region, whereas N r is large in the downwind region. N c is large in the upper layer and is far outside the cell center in R01, which may be due to the presence of a neighboring cell. In contrast, N r is large in the middle levels as well as at the peak altitude in R07. Figure 10a shows the radial profiles of the vertical difference in radiation flux between the surface and cloud-top level D F R , which are averaged around the detected cell center, for the four selected regions. The flux difference increases radially up to around a radius of 0.9 km and then decreases. The magnitude of the radial decrease is more significant in downwind regions. The peak of the flux difference at a radius of 0.8 km is consistent with the radius-height cross sections of q c and q r , in which the column-integrated value is large at the radius.
The averaged cloud depth decreases in the downwind side (Fig. 10b), which is consistent with the difference in radiation flux averaged in the cloud cells. As indicated by the horizontal LWP fields and radius-height sections, the cloud cells in the upwind side have thicker clouds than those in the downwind side. Since the thickness of the clouds does not largely differ among the subdomains (cf. Fig. 4), the flux difference mainly results from the liquid water content.  This results in the variation of the vertical radiation flux in the streamwise direction as shown later.
The composite analysis shows that the structures of the cloud cells are qualitatively the same across the entire simulation domain, whereas the horizontal cloud field (Fig. 1) suggests that the cell structure changes from closed to open in the domain. This is consistent with Sato et al. (2015b), and the simulated cloud field has an open cell-like structure across the entire domain. We interpret that the horizontal cloud field is the result of cloud cells for which the distance between cells is increasing in the streamwise direction (Ovchinnikov et al. 2013;Sato et al. 2015a).
Satellite observations and numerical simulations in previous studies have suggested that the horizontal cloud field of low-level clouds consists of a number of cloud cells with the same structure (Wood and Hartmann 2006;Comstock et al. 2007). Figure 11 shows radius-height cross sections of the standard deviation of q c and q r . The quantities are an order of magnitude smaller than the ensemble mean (cf. Fig. 8). This is the same for other quantities, such as number densities (figure not shown). In conclusion, the simulated cloud field consists of a number of cloud cells with almost the same structure.
To compare the area-integrated vertical fluxes (cf. Fig. 4) with those in the cloud cells, vertical profiles of the quantities in the cloud cells are shown in Fig.   12. The quantities were averaged (Figs. 4a -d) and integrated (Figs. 4e -h) within the 5-km radius from the center of the detected cloud cell. The vertical profiles are qualitatively the same as those averaged in the subdomains, and the order of magnitudes are also the same. Thus, the orders of total vertical fluxes for the quantities are equal to the total fluxes in the cloudy region around the detected cloud cells, whereas the profiles fluctuated. Furthermore, the differences among the four subdomains are qualitatively the same as those of the area-integrated values (Figs. 4e -h). In other words, the vertical fluxes averaged over the subdomain can be represented by the fluxes averaged in the cloud cells.

Energy budget analysis
We performed an energy budget analysis for each subdomain. The conservation equation for moist enthalpy was first derived from the temperature equation and the conservation equation for water vapor by assuming incompressibility, as follows: Fig. 11. Radius-height sections of standard deviation of the mixing ratio of (top) cloud water and (bottom) rain water for R01, R03, R05, and R07, which are obtained from all samples from the cloud cells. Units: kg kg −1 .
where t is time, x j ( j = 1, 2, 3) = (x, y, z) are the spatial dimensions in a Cartesian coordinate system, u j ( j = 1, 2, 3) = (u, v, w) is velocity in the three directions, T (x, y, z, t ) is temperature, q v (x, y, z, t ) is the water vapor mixing ratio, c p is the specific heat at a constant pressure for dry air, L v is the latent heat, K T (x, y, z, t ) is the thermal diffusion coefficient, K q (x, y, z, t ) is the diffusion coefficient for the water content, Q (x, y, z, t ) is the tendency of the water vapor associated with the phase change, D s (z) is the divergence due to large-scale subsidence (s −1 ), U is the background flow speed (m s −1 ) and is constant across the entire domain, F R (x, y, z, t ) is the vertical radiative flux of long wave radiation, and F D (x, y, z, t ) is the vertical flux of dissipative heating. The tensor notation follows Einstein's summation convention. The lefthand-side (LHS) represents the local tendency and advection terms, and the terms on the right hand side (RHS) indicate diffusion, diabatic heating associated with phase change, large-scale divergence, advection due to the large-scale constant flow, radiative flux divergence, and frictional dissipation flux divergence, respectively. Later on, the dissipative heating term is Vertical profiles of (upper) temperature (K), water vapor mixing ratio (g kg −1 ), cloud water mixing ratio (g kg −1 ), rain water mixing ratio (g kg −1 ), liquid-water potential temperature (K), and (lower) vertical fluxes of temperature (K m s −1 ) and water mixing ratios (g kg −1 m s −1 ). The quantities in the upper panels are averaged and the fluxes are integrated within the 5-km radius from the center of the detected cloud cells. They are averaged from t = 12 h to 16 h. The dashed lines in the panels, (a), (b), and (e), indicate the initial values.
neglected in this study, because it is not included in the present simulation. By assuming K T = K q and combining the equations after c p ´ (5) and L v ´ (6), we obtain the conservation equation for moist enthalpy (k = c p T + L v q v ) as The equation was integrated over a volume of subdomain from the surface to the cloud top height as where V and A are the volume and horizontal area of the subdomain, respectively, F sk is the surface flux of moist enthalpy into the atmosphere,  LS (º D s x 3 ¶ x 3 k) represents the large-scale subsidence term, and  CF (º -U ¶ x 1 k) represents the horizontal advection by constant flow. The first and second terms on the RHS are derived from the diffusion term and large scale divergence term in (7). We assumed that advection terms by local velocity (second term on the LHS of (7)) vanish in the integration due to the periodicity at the lateral boundaries of the subdomain, and we assumed that the enthalpy flux at the surface was F sk and vanishes at the cloud top. The integrated equation (8) indicates that the tendency of the enthalpy in a subdomain is determined by the imbalance among the surface flux, large-scale subsidence, constant flow advection, and radiation. Figure 13 shows the budgets of (8) calculated from the simulation averaged from t = 12 h to 16 h. The surface enthalpy flux is the dominant term that contributes positively to the system. The other terms play a negative role, and their sum corresponds approximately to the surface enthalpy flux. The radiation flux outside the cell region is nearly constant in the x direction, while the magnitude of the radiation flux in the cloud cells increases. The magnitude of large-scale subsidence is not large, but it is not negligible either, and increases toward the downwind side. The magnitude of the advection associated with the constant background velocity is the largest among the negative terms, and this is caused by a combination of the background flow (~ 5 m s −1 ) and the spatial enthalpy gradient. The surface flux that is large in the downwind region increases enthalpy in the air and enhances the spatial gradient of enthalpy and the horizontal advection term, while it does not largely change in the x direction.

Development of a cloud fraction model
Now we assume that the cloud field in a subdomain consists of a finite number of cloud cells with the same structure, and also that the radiation flux in the energy equation is considered separately in the cloudcell region and other areas. By assuming that the cloud cell has a radius of r c , (8) can be written as where A is the area of the subdomain, A c is the total area of cloud cells defined by Nò 0 rc 2πrdr, N is the number of cloud cells in a subdomain, the tilde denotes the average inside of the cloud cell, the hat denotes the average outside of the cloud cell, and D is the vertical difference between the surface and the cloud top level. This equation can be written as á F sk ñA + á  LS ñ AH + á  CF ñ AH where H is the depth of the cloud layer and the angle bracket indicates the domain average. We have assumed a quasi-steady state (á ¶ t k ñ ® 0). This results in an equation for the cloud fraction as follows:Ĉ Thus, the cloud fraction is determined from a balance between the surface enthalpy flux with large-scale subsidence and constant flow advection, and radiation fluxes in cloud cells and other areas. As the vertical difference in the radiation flux represents a flux divergence that tends to remove energy from the domain, they are negative and the minus sign on the RHS of (11) is reasonable unless the negative contribution of large-scale subsidence and the advection terms are greater than the energy input from the surface enthalpy flux. This would happen only rarely because the surface flux is on the order of 100 W m −2 , whereas the large-scale subsidence and advection terms are less than 100 W m −2 , or rather, an order of magnitude smaller than the surface flux in the focusing regions (cf. Fig. 13). The model implies that the cloud fraction increases with a large surface flux, large radiation flux, and (negatively) small effects of background flow motion. This suggests that low-level clouds cannot be present (because the cloud fraction would become less than 0) if the large-scale forcing is stronger than the surface flux; i.e., when the energy removal by largescale forcing exceeds the energy input to the system from the ocean. Figure 14 shows the cloud fractions (temporally averaged from t = 12 h to 16 h) that were estimated from our simulation and diagnosed using (11) for each subdomain. The simulated cloud fraction monotonically decreases in the streamwise direction from 0.99 in R01 to 0.66 in R09 as shown in Figs. 1 and 3.
The diagnosed cloud fraction decreases in the streamwise direction, whereas it slightly overestimates the simulated fraction especially in the downwind regions. The standard deviation of CF, which is obtained in each subdomain from t = 12 h to 16 h, is large in the upwind regions. It is implied that the standard deviation of radiation fluxes are small in the upwind regions as the cloud field appears not to largely change in time (cf. Fig. 3). Since the radiation fluxes are in the denominator of (11) and there is no temporal variation in the surface flux in the simulation, which has the largest values among the terms in (11), it is suggested that the standard deviation of CF diagnosed by (11) is large in the upwind regions.
The relative contribution of the terms of (11) is listed in Table 1. On average, the surface flux term is the source of energy and the other terms remove the energy from the volume (cf. Fig. 13). In order to examine the sensitivity of the cloud fraction to individual terms in (11), we calculated the cloud fraction by artificially increasing the magnitude of each term by 10 % from the spatiotemporal average that was used to calculate the mean value in Fig. 14, while the others were fixed. Figure 15 shows the diagnosed cloud fraction by (11) with the magnitude of each term increased individually by 10 %. The cloud fraction is most sensitive to the surface flux, which is followed by the horizontal advection term by constant flow and the longwave radiation term in cloud cells. Although the magnitudes of change are large, it should be noted that the sensitivity has been tested by artificially increasing a single term by 10 %, while the other terms are fixed. This never happens in nature because a change in one term more or less affects other terms.
The results are summarized in Table 1, which lists the rate of change in the calculated cloud fraction. It indicates that a term with a larger mean value has a larger sensitivity as expected; the largest sensitivity of cloud fraction is found in the surface flux term, which Fig. 14. Cloud fraction estimated from the simulation (cf. Fig. 3) (black line) and diagnosed using Eq. (11) (gray line). The data are temporally averaged from t = 12 h to 16 h. The gray hatched area stands for the ± 1 standard deviation from the mean, which are estimated using quantities spatially averaged in each subdomain from t = 12 h to 16 h.
has the largest mean value in (11). The large contribution of horizontal advection by constant flow implies the importance of representation of the wind direction as well as the wind speed in global models, while the variance of wind direction is large in models  in the Coupled Model Intercomparison Project (CMIP) (Taylor et al. 2012).
Another balance equation can be obtained by using moist conservative variables such as the liquid water potential temperature θ l , and the concept of the present study can be applied. In this study, we consider the moist enthalpy, as it is easily obtained from T and q c , which are the prognostic variables in the numerical model used for the present simulation.

Discussion
One of the applications of the model developed here to diagnose the fraction of low-level clouds is to use it as a parameterization when estimating the radiation budget of global models. The model is valid when a sufficient number of cloud cells is included in a domain that corresponds to a grid point in the global model. Specifically, the present model was developed for an area of 30 ´ 28 km 2 , which is smaller than the grid spacing of conventional climate models.
The diagnostic equation has been simplified using some assumptions. Using quantities predicted or diagnosed in numerical models, terms in the numerator on the RHS of (11) (i.e., surface enthalpy flux, large-scale subsidence, and large-scale constant flow) can be calculated. Nevertheless, radiation fluxes in the cloud cell and other areas in the denominator are not explicitly estimated by climate simulations, as parameterizations are needed to calculate the fluxes from grid-scale quantities. Thus, some additional assumptions are necessary to estimate the fluxes. Once the differences in radiation fluxes are approximated using the gridscale quantities from the climate models, the cause of bias in the cloud fraction in the numerical models can be detected by comparing the estimated terms in (11) with the outputs of climate simulations.
A possible approach to estimating the vertical difference in radiation fluxes, especially the radiation flux in the cloud cell region, would be to use LWP. The present simulation uses a parameterization for radiation flux based on the vertically integrated q l , suggesting that LWP could be a useful quantity. For example, Chung et al. (2012) proposed a method to evaluate the radiation flux term with the usage of the probability density function (pdf) of LWP. As they applied a gamma function for the pdf of LWP, the radiation flux is a function of the mean, the homogeneity, and the standard deviation of LWP. Their formulation could be applied to the present equation for the radiation flux and hence the cloud fraction. However, accurate estimation of LWP in the lower layer in climate models is also a difficult issue.
It should also be noted that the parameterization of the radiation flux (Stevens et al. 2005a) used in the present simulation is designed for nocturnal longwave radiation. Hence, a more realistic radiation scheme, such as one including short-wave radiation, would be required to apply the model to more realistic cases.  The budget equation (10) would possibly be more accurate, if the diffusion term or eddy vertical advection term at the top of boundary layer are considered, which have been neglected to derive (10). As introduced above, cloud top entrainment and shallow convection would also play important roles in determining structures of low-level clouds and these processes would be incorporated in the two terms. Further studies are needed to formulate the effects of processes at the cloud top as an area-integrated flux.
Additionally, considering the tendency term á ¶ t kñ or the advection term may also help to increase the accuracy of the CF equation (11), whereas they have been neglected as we assumed the steady state and periodicity in each subdomain. The CF shown in Fig. 14 was estimated by quantities that are averaged in a subdomain, not in time, and then the average and variance of CF were calculated. Hence, the assumption of a quasi-steady state would not work and it is desired in a future study to discuss the statistics of diagnosed CF by using different data sets, each of which individually with different quasi-steady state.
It is worth discussing the physical reason why the diagnosed CF by (11), which is derived from an energy balance equation, decreases in the streamwise direction. Let us consider a simple system in which each terms in (11) change linearly in the streamwise direction. In this case, (11) can be re-written as where α is the proportional coefficient for surface flux in the x direction, β is the proportional coefficient for all the other terms in the numerator of (11), γ is the coefficient for the terms in the denominator, and B 1 and B 2 are the values of the terms in the numerator and denominator at the upwind edge. Taking a spatial derivative in the x direction yields The condition to decrease CF in the x direction ( ¶ x C F < 0) can be written as provided that γ > 0 and γx + B 2 > 0, which is implied by Fig. 13. The RHS of the inequality corresponds to CF (12), which is approximately 1 or less. Thus, when the ratio of coefficients in the LHS is roughly less than 1, the CF tends to decrease in the x direction. The result of the simulation (cf. Fig. 13) shows that α is positive and the largest contributor, β is small and negative, and γ is positive and large, which possibly results in the decreasing trend of CF. A large α tends to increase the CF as implied in (11), whereas large β and γ result in a decreasing trend of CF. In other words, gradual increase in surface flux and rapid decreases in large-scale effects and radiation fluxes are favorable for decrease in CF. An increase in surface flux would enhance the LWP, which increases radiation flux in cloud cells per unit area DF  R . Since the water vapor is larger in the downwind side, a larger amount of latent heating appears to be released once condensation occurs. It is possible under an environment in which net income fluxsurface flux together with the other large-scale effects that are in the numerator of the LHS of (14) and tend to decrease the increasing rate of surface fluxincreases in the streamwise direction. Larger latent heating more likely produces a positive feedback between the heating and the vertical motion, which results in a larger amount of condensed water in the vertical direction and also in a smaller area (Bjerkness 1938). It increases DF  R in the streamwise direction and hence the increasing rate of radiation flux may exceed that of the net income flux. Thus, condition (14) can possibly be satisfied, indicating a decrease in diagnosed CF in the streamwise direction.

Conclusions
In this paper, we have proposed an energy balance model to diagnose the fraction of low-level clouds using a conservation equation of moist enthalpy and by assuming that the horizontal field of low-level clouds consists of a number of cloud cells with the same structure. The derived equation indicates that the cloud fraction can be represented as the ratio of the sum of the surface enthalpy flux, large-scale subsidence, and large-scale constant flow to the radiation fluxes in the cloud cells and other regions.
Energy budget analysis was performed on the results of a simulation covering a wide area and with fine spatial resolution. The surface flux played a dominant role in transferring energy into the domain, while the other terms (radiation, large-scale subsidence, and advection by large-scale motion) made a negative contribution, and their magnitudes were small. Detecting cloud cells in the simulation showed that the structure of cloud cells is qualitatively the same in all subdomains. Furthermore, the standard deviation of the quantities of cloud cells was smaller than the mean values, and the total fluxes in each subdomain were of the same order of magnitude as the total fluxes transported by all of the detected cloud cells. This indicates that the cloud field of low-level clouds consists of a number of cloud cells with the similar structure, which is consistent with previous studies.
We developed the cloud fraction model based on the following assumptions: the low-level cloud field consists of a finite number of cloud cells, each of which has the same structure, and the radiation fluxes are separately considered in cloud cell field and other regions. In theory, the model is applicable to low-level clouds with both closed and open structures as long as the cloud field consists of a number of cloud cells with almost the same structure. The model was tested to diagnose the cloud fraction of the simulation, and it was able to quantitatively capture the cloud fraction reasonably well. Our results indicate a potential use of this model for parameterization in climate models, although a tuning process and the closing of the equation would be required. Further verifications are needed to apply the developed parameterization to global simulations.