444 MONTHLY WEATHER REVIEW VOLUME 139 Interpolation of LES Cloud Surfaces for Use in Direct Calculations of Entrainment and Detrainment JORDAN T. DAWE AND PHILIP H. AUSTIN Department of Earth and Ocean Sciences, University of British Columbia, Vancouver, British Columbia, Canada (Manuscript received 22 April 2010, in final form 4 August 2010) ABSTRACT Direct calculations of the entrainment and detrainment of air into and out of clouds require knowledge of the relative velocity difference between the air and the cloud surface. However, a discrete numerical model grid forces the distance moved by a cloud surface over a time step to be either zero or the width of a model grid cell. Here a method for the subgrid interpolation of a cloud surface on a discrete numerical model grid is presented. This method is used to calculate entrainment and detrainment rates for a large-eddy simulation (LES) model, which are compared with rates calculated via the direct flux method of Romps. The comparison shows good agreement between the two methods as long as the model clouds are well resolved by the model grid spacing. This limitation of this technique is offset by the ability to resolve fluxes on much finer temporal and spatial scales, making it suitable for calculating entrainment and detrainment profiles for individual clouds. 1. Introduction The largest uncertainties in climate sensitivity estimates from general circulation model (GCM) simulations come from the subgrid-scale parameterization of low clouds (Colman 2003; Bony and Dufresne 2005; Webb et al. 2006). Specifically, stratocumulus clouds and the transition regime from stratocumulus to trade cumulus are the dominant source of variance between models in the estimation of the cloud radiative response to changing climate (Williams and Webb 2009). Proper simulation of the subgrid-scale effect of cumulus clouds in GCMs requires understanding the rates at which air is entrained into and detrained from the clouds. Cloud entrainment and detrainment rates exert influences on profiles of cloud properties, the height of the cloud tops, the amount of heat and moisture the clouds transport upward, and the heights at which the clouds deposit that heat and moisture. They also have effects on the vertical transport of aerosols out of the boundary layer and the rate at which chemical reactions can occur in those aerosols (Barahona and Nenes 2007; Andrejczuk et al. 2008). Corresponding author address: Jordan T. Dawe, Department of Earth and Ocean Sciences, University of British Columbia, 6339 Stores Rd., Vancouver, BC V6T 1Z4, Canada. E-mail: jdawe@eos.ubc.ca DOI: 10.1175/2010MWR3473.1 Ó 2011 American Meteorological Society Large-eddy simulation (LES) is an important tool used in the study of cloud entrainment and detrainment. LES models apply grid resolutions on the order of 10–100 m, well within the domain of the Kolmagorov 25/ 3 turbulence spectrum. This allows for a relatively simple turbulence model that captures the important statistics of the subgrid-scale eddy fluxes and thus, an accurate representation of the atmospheric physics in a domain ’10 km2, which is large enough to simulate a field of shallow clouds. LES can be assessed against results taken from field surveys such as the Barbados Oceanographic and Meteorological Experiment (BOMEX; Holland and Rasmusson 1973) or the Atmospheric Radiation Measurement Program (ARM; Brown et al. 2002), and such comparisons show good agreement between LES and data. Several recent studies have looked at the life cycle of individual clouds taken from LESs, trying to break the cloud field into its component parts (Zhao and Austin 2005a,b; Heus et al. 2009). Estimates of entrainment and detrainment rates for individual clouds would be quite useful in these studies, but are difficult to achieve. Entrainment and detrainment rates for nonprecipitating clouds are typically calculated in LES by recording budgets of bulk conserved tracer variables, such as the total humidity or the liquid water moist static energy, and inferring the amount of fluid exchange between cloud and clear air that is needed to explain the rate at which that tracer is being vertically advected within the cloud field FEBRUARY 2011 (Siebesma and Cuijpers 1995). These budgets typically assume the clouds and the cloud environment are horizontally homogeneous slabs; this is a much less accurate assumption on the level of an individual cloud, and this variability makes bulk tracer calculations on this scale problematic. Alternatively, entrainment and detrainment could simply be calculated directly from the LES velocity, pressure, humidity, and temperature fields. Siebesma (1998) defines entrainment and detrainment as E 5À D5 445 DAWE AND AUSTIN 1 A 1 A þ ^Á(uÀui ),0 n þ ^ Á(uÀui ).0 n r^ n Á (u À ui ) dl, r^ n Á (u À ui ) dl, (1) (2) where E and D are the entrainment and detrainment rates (in kg m23 s21), r is the density of air (in kg m23 s21), u is the velocity of the air (in m s21), ui is the velocity of the cloud surface (in m s21), A is the area of the cloud (in m), ^ is a unit vector directed out the cloud surface, and the n path integral is taken around the cloud surface at a constant vertical level. However, the accuracy of this method applied to LES suffers from the need to calculate the velocity of the air relative to the cloud surface. In reality these velocities are very nearly identical, ;1–2 m s21, but the discrete nature of the LES model grid forces the modeled surface velocity to be either 0 m s21 or Dx/Dt ’ 15– 30 m s21, where Dx is the model grid spacing and Dt is the model time step. The surface of the cloud only moves when a grid cell’s humidity reaches saturation, and when it does, an entire gridcell worth of fluid leaves or enters the cloud. Romps (2010) recently overcame this problem by averaging the cloud surface fluxes over the time needed for an entire grid cell to entrain or detrain. This is done by summing the ‘‘activity source’’ at the cloud boundary, › (rA) 1 $ Á (ruA), ›t (3) over the period in which the cloud surface is adjacent to a grid cell. Here A is the ‘‘activity’’ of the air, a scalar field that is 1 inside the cloud and 0 outside. The activity can be easily redefined to determine entrainment and detrainment for any material surface in a numerical model; Romps used a definition similar to Siebesma and Cuijpers (1995)’s cloud core criteria (fluid with condensed liquid water, positive buoyancy, and upward velocity) in his calculations. Once a grid cell is no longer adjacent to the boundary between active and inactive air, the sign of the summed activity source is evaluated. If the summed activity source is positive, more mass has been entering the region of active air than has been leaving, and the summed activity source is considered to be entrainment. Conversely, a negative sum is considered to be detrainment. (In this paper, we shall refer to this method of summing the activity source to determine entrainment and detrainment as the ‘‘Romps method.’’) Romps found entrainment and detrainment values resulting from this method were approximately twice as large as those from previous bulk tracer calculations, suggesting bulk tracer calculations contain significant biases. However, fine temporal and spatial resolution of entrainment and detrainment variability is difficult to achieve with this method because of the long temporal averages needed. Here we present a method for calculation of the cloud entrainment and detrainment rates that relies on interpolation of the subgrid location of the cloud surface. This method can be used to produce accurate estimates of the cloud entrainment and detrainment rates for individual LES clouds. In section 2 we describe this method, in section 3 we describe the model configuration we used to test this method, in section 4 we compare this calculation with entrainment and detrainment rates calculated using Romps (2010) direct flux calculation, and in section 5 we discuss our results and present our conclusions. 2. Method Here we derive a method for the calculation of entrainment and detrainment fluxes through the surface of an arbitrary volume in a numerical model. In our work this is used to calculate the mass exchange between cloud and environment, but this method could be used to calculate mass flux into and out of any volume of interest. To maintain generality we do not yet specify how the location of this cloud surface is determined. Consider a numerical model grid cell containing a cloud surface with normal vector C, where C points outward from the surface and has units of meters squared (Fig. 1). This surface, combined with W, the portion of the grid cell walls that lie within the cloud (in m2), encloses a cloud volume V that has units of meters cubed. Siebesma (1998) gives the net entrainment and detrainment over the cloud surface to be E À D5 ð C r(ui À u) Á dC, (4) where r is the air density (in kg m23), u is the velocity of the air (in m s21), and ui is the velocity of the cloud interface (in m s21). Calculating this integral requires knowledge of the velocity field over the surface of C and the time evolution of C, neither of which is easily calculated in a numerical model. Instead, we seek a simplified but equivalent calculation. 446 MONTHLY WEATHER REVIEW VOLUME 139 E À D5 r d dt ð dV À V (t ) ð ru Á dC. (7) C Next we apply the divergence theorem to simplify the flux integral through C: ð $ Á (ru) dV 5 V ð ru Á dC 1 C ð ru Á dW. (8) W As a result of the mass conservation, $ Á (ru) 5 0, which implies ð ð ru Á dC 5 À ru Á dW. C FIG. 1. Gridcell schematic showing geometry for the derivation of direct entrainment and detrainment fluxes using a subgrid-scale interpolation of the cloud surface location. Vector C is normal to the arbitrary cloud surface, shown by the thick line. Vector W is normal to the portion of the gridcell walls occupied by cloud, shown by the dotted lines; Wx is used to denote the x direction grid cell surfaces; and Wy denotes the y direction surfaces. These surfaces enclose a volume V, shown by the shaded area. The u and y vectors represent the locations of the horizontal velocities in the x and y directions, respectively, for an Arakawa C-grid. To calculate the velocity of the cloud surface, we make use of the Leibnitz theorem: ð ð ð d ›r dV 1 r dV 5 rui Á dC dt V (t) V (t) ›t C(t) ð rui Á dW. (5) 1 W (t ) Since the walls of the grid cell do not move, ui is 0 over W. If we also assume ›r/›t ’ 0, we get r d dt ð ð dV 5 V (t ) C (t ) rui Á dC. (6) We then combine (4) and (6) to give Substituting this into (7) results in E À D5 r dV 1 dt ð ru Á dW. (10) W Following Romps (2010), if (10) is positive, the result of this calculation is assumed to all be E, and if negative, D is ð dV 1 ru Á dW , E 5 max 0, r dt W (11) ð dV À ru Á dW . D 5 max 0, Àr dt W (12) Therefore, we can find the entrainment and detrainment by calculating the rate of change of the cloud volume inside the grid cell and the mass flux through the cloudy portion of the gridcell walls. We assume the mass flux is constant over the gridcell wall; this makes calculating mass fluxes trivial for a model using an Arakawa C-grid, while interpolation can give fluxes for other grid configurations. Expanding the terms in (10) for an Arakawa C-grid-type model with a horizontally uniform density gives dV (i, j,k) 1 r(k) u(i11, j,k) W x(i11, j,k) À r(k) u(i, j,k) W x(i, j,k) 1 r(k) y (i, j11,k) W y(i, j11,k) dt À r(k) y (i, j,k) W y(i, j,k) 1 1/2[r(k11) 1 r(k) ]w(i, j,k11) W z(i, j,k11) À 1/2[r(k) 1 r(kÀ1) ]w(i, j,k) W z(i, j,k) . E À D 5 r(k) (9) W There are a multitude of interpolation schemes that can be used to determine the cloud volume and cloudy grid wall area in a numerical model. The simplest would be to assume that cloudy grid cells are completely filled with cloud, and otherwise grid cells contain no cloud at all. We refer to this as the ‘‘no interpolation’’ method. Since the cloud surface only moves when a whole grid (13) cell undergoes condensation or evaporation this will result in poor estimates of dV/dt, C, and W, which will result in an overestimate of the entrainment and detrainment. To reduce this overestimate, we interpolate the location of the cloud surface at each time step. Next we describe the interpolation scheme we use for this purpose. FEBRUARY 2011 447 DAWE AND AUSTIN FIG. 2. Schematic representation of the tetrahedral interpolation method. (left) The subdivision of the model grid cell into 48 tetrahedrons. (middle) The three cases for interpolation of the cloud surface within a subtetrahedron between points at the four vertices of the tetrahedron; all other tetrahedron interpolation cases are variations on these three. Black dots represent cloudy tetrahedron points, white dots represent clear tetrahedron points, and gray dots represent the interpolated cloud position between cloudy–clear tetrahedron point pairs. (right) A horizontal slice through the model’s Arakawa C-grid. Rightward-pointing triangles represent u-velocity points, upward-pointing triangles represent y-velocity points, circles represent bulk tracer quantities such as temperature or humidity, and x represents the location of the cloud surface interpolation points. Dotted lines show the boundaries of the tetrahedral subdivision of the grid, and shading represents cloud volume. Cloud surface interpolation Several standard techniques exist for subgrid isosurface interpolation in the field of computer visualization, such as the Marching Cubes algorithm (Lorensen and Cline 1987). In the spirit of these techniques, we have implemented a scheme that relies on subdividing the grid cells into regular subcells, finding the location of the cloud surface within these subcells, and then reassembling the cells to construct the location of the cloud surface. Our interpolation scheme, which we refer to as ‘‘tetrahedral interpolation,’’ splits the cubic grid cell into six pyramids, then splits each pyramid into eight tetrahedrons (Fig. 2, left panel). This results in 48 tetrahedrons, each composed of 4 vertices located at the gridcell center, the center of a gridcell wall, the center of a gridcell edge, and a gridcell corner. Generally, model scalar quantities will be defined on only one of these vertices, and so we must interpolate the surrounding grid to determine values at the other three points. On an Arakawa C-grid, humidity is defined at the grid centers, and so we must interpolate to find the values at three of the tetrahedron vertices. To find values at the center of the gridcell wall, we average the two values on opposite sides of the wall; at the gridcell edge, we average the four values around the edge; and at the corners, we average the values of the eight cells surrounding each corner. The cloud surface itself is defined in a similar way to the Marching Cubes algorithm (Lorensen and Cline 1987). Each tetrahedron has four vertices, which can either be inside or outside of the cloud. This results in 16 possible combinations of vertex values that must be considered. Many of these cases share symmetries, reducing the number of independent cases to three classes. If only one of the tetrahedron vertices is cloudy (or conversely, only one is clear) than that corner of the tetrahedron is cut by a triangle representing the cloud surface (Fig. 2, center panel, upper left); this accounts for eight of the cases. If two vertices are cloudy and two vertices are clear, this results in the tetrahedron being cut by two triangles that share a common edge (Fig. 2, middle panel, top right); this accounts for six of the cases. Finally, if the vertices of a tetrahedron are all cloudy or clear, the surface does not pass through the tetrahedron (Fig. 2, middle panel, bottom); this accounts for the remaining two cases. The locations of the vertices of the triangles that cut the tetrahedron are found by linear interpolation between cloudy and clear tetrahedron vertices. Repeated application of this algorithm to each grid cell results in a continuous surface that approximates the subgrid location of the cloud surface (Fig. 2, right panel). Once the geometry of the case is determined, the area of the cloudy portion of the model gridcell wall can be calculated by dividing the cloudy area into triangles and summing the triangle areas using the well known crossproduct area formula: A5 j(b À c) 3 (a À c)j , 2 (14) where a, b, and c are the position vectors of the vertices of the triangles. Similarly, the volume of the cut tetrahedron is calculated by subdividing it into smaller tetrahedrons, and summing the volumes of the subtetrahedrons using the triple-product volume formula (Harris and Stocker 1998): 448 MONTHLY WEATHER REVIEW VOLUME 139 FIG. 3. Comparison of mean model profiles of core (a) entrainment and (b) detrainment calculated using the tetrahedral interpolation method (black line), the Romps method with upwind advection (dotted line), and the Romps method with MPDATA advection (gray line). (d),(e) As in (a),(b), but for the fractional entrainment and detrainment rates. V5 j(a À d) Á [(b À d) 3 (c À d)]j , 6 (15) where a, b, c, and d are the position vectors of the vertices of the subtetrahedrons. This results in all the information needed to measure entrainment and detrainment via (10). 3. Model description We have implemented both our tetrahedral surface interpolation scheme and the time averaging scheme of Romps (2010) in the System for Atmospheric Modeling (SAM; Khairoutdinov and Randall 2003), allowing the model to calculate entrainment and detrainment via these two methods simultaneously. In this paper we perform these calculations for the cloud core, which we define following Siebesma and Cuijpers (1995) as fluid having condensed liquid water, upward velocity, and positive buoyancy relative to the horizontal mean. The cloud core surface then is located where the total humidity equals the saturated humidity, Ty is equal to the horizontal mean of Ty, and the vertical velocity is zero. In the tetrahedral method, finding the location where qt equals qsat (T, p) is complicated by the nonlinearity of qsat (T, p). For simplicity, we define a variable qdiff 5 qt 2 qsat (T, p), where qt is the total specific water content and qsat (T, p) is the saturated specific humidity, all in FEBRUARY 2011 DAWE AND AUSTIN kilograms of water per kilograms of moist air. Here qdiff can be thought of as the moisture excess or deficiency relative to the saturated humidity and qdiff 5 0 is the location where the total humidity equals the saturated humidity. Similarly, we calculate DT y 5 T y À T y , where the bar represents the horizontal mean, and Ty is the virtual temperature in kelvin, to represent buoyancy excess or deficiency. Vertical velocity w requires no new variable to be calculated, but as SAM is an Arakawa C-grid model, it is defined at the center of the top and bottom of the grid cell instead of at its center. Therefore, we first interpolate w to the center of the grid cell for consistency with location of the specific humidity and virtual temperature. To find the location of the cloud core between a point that is in the core and a point that is outside the core, the locations between these points where qdiff, DTy, and w are 0 are calculated, and the value that is closest to the cloud core point is selected as the start of the region where all the cloud core criteria are satisfied. Since the Romps scheme relies on calculating the advection of active (i.e., cloud core) air, there are actually several ways this calculation can be performed. Romps (2010) used a simple first-order upwind advection scheme; we have implemented both this scheme and one based upon SAM’s second-order Multidimensional Positive Definite Advection Transport Algorithm (MPDATA) routine (Smolarkiewicz and Grabowski 1990). We have run these schemes in a standard Global Energy and Water Cycle Experiment Cloud System Study BOMEX LES (Holland and Rasmusson 1973; Siebesma et al. 2003). The model was run with a domain extent of 6.4 km 3 6.4 km in the horizontal and 3.2 km in the vertical. Model cloud area, vertical velocity, and cloud core entrainment diagnosed from bulk conserved tracer budgets agree well with the results presented in Siebesma et al. (2003) (not shown). To examine the resolution dependence of our scheme, we ran three models: one at 25-m grid spacing in all directions with a 1.5-s time step, one at 50-m grid spacing with a 3-s time step, and one at 100-m grid spacing with a 6-s time step. For most of the results we present here we rely on the 25-m resolution model. The models were each run for 6 h, and the first 3 h of the simulation were discarded as the model was still approaching the steady state. The 15-min averages were output for the terms of each of our calculations. Running the 25-m resolution model with the tetrahedral surface interpolation scheme resulted in a 13% model slow down. 4. Comparison with Romps Here we compare our tetrahedral interpolation method with the Romps method calculated using upwind advection 449 FIG. 4. Average profiles of entrainment minus detrainment calculated using the continuity equation (dotted line), direct fluxes with no surface interpolation (gray line), the tetrahedral interpolation method (thick black line), and the Romps method using upwind advection (thin black line). and MPDATA advection. Comparison of the average values produced by the tetrahedral interpolation method and the Romps method with MPDATA advection over 3 h of model time (Fig. 3) shows good agreement in the vertical profiles of both E and D; however, the magnitude of the tetrahedral entrainment and detrainment values are about 15% higher than the Romps MPDATA value. The Romps method with upwind advection, however, gives results that are about 50% larger than the tetrahedral method. This difference is likely the result of larger numerical errors in the first-order upwind scheme compared with the MPDATA scheme, and is therefore an estimate of the influence of the numerical method used on the values of E and D. In either case, the tetrahedral scheme is vastly superior to using no interpolation, which results in entrainment and detrainment values four times larger than the Romps MPDATA value (not shown). There are similar levels of agreement between 5 E/M and d 5 D/M, the fraction entrainment and detrainment rates, where M is the total vertical mass flux within the cloud core (in kg m22 s21) and and d have units of per meters. Since the tetrahedral method interpolation redefines the extent of the cloud core, we calculate separate values of M for the Romps method and the tetrahedral method. This results in slightly smaller M values for the tetrahedral calculation. Both and d show large deviations between the tetrahedral and Romps calculations 450 MONTHLY WEATHER REVIEW VOLUME 139 FIG. 5. Profiles of the correlation between the Romps method and the tetrahedral interpolation method for (a) and (b) d. Black lines show the correlation using upwind advection in the Romps method and gray lines show the correlation using MPDATA advection. The dotted line shows the 99% confidence level for a significant correlation. near the cloud base and top. These deviations likely result from the small cloud fractions at these heights, which both makes the statistical sampling less robust and allows small differences between E and D to be magnified by small values of M. a. Agreement with the continuity equation While E and D are important for calculating changes in cloud properties such as temperature and humidity, the vertical profile of cloud fraction is controlled by their difference, E 2 D. For comparison with these direct flux methods, we also calculate E 2 D from the continuity equation for a turbulent plume: r ›a ›M 1 5 E À D, ›t ›z (16) where r is the density of air (in kg m23), a is the cloud fraction, and M is the fractional vertical mass flux (in kg m22 s21). To satisfy continuity, the difference between the amount of fluid entrained and detrained by the clouds at a given height must equal the vertical gradient in cloud mass flux plus the rate in change of the cloud area. We calculate the vertical gradient of mass flux in (16) by interpolating vertical velocities from the edges to the center of each grid cell, then averaging total mass flux within the clouds using these interpolated velocities. The 15-min-average values of mass flux were output and the vertical derivative taken via centered differencing. Also, the 15-min-average values of ›a/›t are calculated within the model. Comparison of the tetrahedral values of E 2 D that we have calculated with the Romps and no-interpolation methods shows reasonable agreement (Fig. 4). Of these three calculations, the no-interpolation method is the most accurate measurement of E 2 D, as we calculate its values directly from the model’s velocity field without any modification of the cloud surface. This is confirmed by the close agreement between the no-interpolation method and the continuity equation. Our direct flux calculations also agree with the continuity equation fairly well (Fig. 4), but diverge between cloud base and 1-km height. The tetrahedral interpolation method is biased low by the method’s redefinition of the cloud volume, while the FEBRUARY 2011 DAWE AND AUSTIN 451 FIG. 6. Time variability of (a),(c) entrainment and (b),(d) detrainment horizontally averaged over the model domain at a height of 1 km calculated via (a),(b) the tetrahedral interpolation method and (c),(d) the Romps method using upwind advection. Romps method’s averaging redefines the time at which entrainment and detrainment events take place. Both these schemes result in a less negative E 2 D value than the nointerpolation scheme, with the tetrahedral interpolation having slightly more bias. Finally, we note that the E 2 D, which results from the Romps MPDATA calculation (not shown), is essentially identical to the Romps upwind calculation. b. Correlation between Romps and tetrahedral methods So far we have shown that our tetrahedral interpolation produces entrainment and detrainment values that agree reasonably well with those produced by the Romps method. However, this does not ensure that the same variability is displayed by the two calculations. To analyze the variability, we take the correlation of the 15-minaveraged and d values calculated via the two methods over the 3-h model run at each height to generate correlation profiles. We do these calculations for both Romps using upwind advection and Romps with MPDATA advection. We use 5 E/M and d 5 D/M to analyze the variability as E and D are both strongly correlated with cloud volume; more cloud volume with the same entrainment velocity results in larger entrainment. Dividing by the mass flux removes this area dependence. Separate M values are calculated for the tetrahedral and Romps methods to take into account the volume redefinition that occurs in the tetrahedral method. Heights at which the model does not have clouds for the entire 3-h period are excluded from the calculation. Correlations between the Romps method, using both upwind and MPDATA advection, and the tetrahedral interpolation method for both and d are significant at the 95% confidence level over most of the cloud layer (Fig. 5). The d shows nearly perfect correlation, while is weaker but still significant. Values near the cloud base and top show no correlation, which we would expect 452 MONTHLY WEATHER REVIEW VOLUME 139 FIG. 7. Vertical cross section of a model cloud showing 1-min-averaged (a),(c) entrainment and (b),(d) detrainment calculated using the (a),(b) tetrahedral interpolation method and (c),(d) the Romps method using upwind advection. Grayscale indicates amount of entrainment or detrainment and lines show the region where condensed liquid water exists. given the poor agreement between the average values of and d between the two methods in these regions (Figs. 3c,d). We take these results to show general good agreement between our tetrahedral interpolation scheme and the time averaging scheme of Romps over relatively large time and space averages. c. Behavior over different time and space scales The Romps method relies on averaging over sufficiently long time scales for several grid cells to complete an entrainment or detrainment cycle. The tetrahedral interpolation scheme, on the other hand, has no such limitation. For example, the tetrahedral interpolation values horizontally averaged over the whole model domain at a given height show little change in variability between instantaneous values, 1-min-average values, and 5-minaverage values (Fig. 6). The Romps method, on the other hand, shows jumps on the order of 50% of the mean from time step to time step. Over 1-min averages, the Romps method has variability similar to that shown by the tetrahedral interpolation method without any averaging, and FEBRUARY 2011 DAWE AND AUSTIN 453 FIG. 8. Mean (top) entrainment and (bottom) detrainment calculated at (left) 25-, (middle) 50-, and (right) 100-m grid spacing using the Romps method with upwind advection (gray line) and the tetrahedral interpolation method (black line). does not settle down to a stable value until 5-min averages are taken. When the spatial variability of the two schemes is examined, the Romps method shows sparse measurement coverage, even when 1-min-average values are taken (Fig. 7). Tetrahedral interpolation values show a relatively smooth spatial distribution at the scale of an individual cloud, and tend to be smaller than the Romps method values. This is the primary advantage of the tetrahedral interpolation method: it allows instantaneous measurements of entrainment and detrainment rates at the individual cloud level. d. Dependence on model resolution Both the Romps and tetrahedral methods display relatively strong dependence on model resolution, with E and D decreasing as grid resolution becomes coarser (Fig. 8). Furthermore, while E and D calculated via the two direct methods show good agreement at 25-m resolution, the tetrahedral interpolation method systematically underestimates the values calculated by the Romps method as grid spacing becomes coarser; at 100-m resolution, the difference in the values exceeds 50%. Since poor estimation of motion of the cloud surface relative to the air should on average result in an overprediction of E and D, a different error must be dominating the tetrahedral calculation. This error is likely related to a systematic underestimate of the cloud volume as calculated by the tetrahedral interpolation method when estimating the volume of a cloud that occupies a single model grid cell. Consider such an isolated cloud cell, surrounded by clear grid cells in the 2D plane (Fig. 9). For simplicity, instead of considering the full cloud core calculation, let us consider simply the condensed water. Assume the cloud cell has a qdiff value of 1 g kg21, and the clear cells, 21 g kg21. The Romps calculation, in agreement with the model assumptions, will treat this entire grid cell as containing cloud. The tetrahedral calculation, on the other hand, will underestimate the cloud volume. At the cell edge points, the qdiff value will be averaged between adjacent cells, giving qdiff 5 0 and placing the cloud surface at the gridcell wall, exactly as expected. At the cell corners, however, the four surrounding qdiff values average to give qdiff 5 20.5 g kg21, and interpolation of this value between the corner and 454 MONTHLY WEATHER REVIEW FIG. 9. Schematic showing bias in the tetrahedral interpolation method, which leads to underestimation of cloud volume. Rightwardpointing triangles represent u-velocity points, upward-pointing triangles represent y-velocity points, and circles represent bulk tracer quantity points. Black circles indicate cloudy points with total specific water 1 g kg21 above the saturated specific humidity and white circles indicate clear points with total specific water 1 g kg21 below the saturated specific humidity. The gray area indicates the resulting tetrahedral interpolation cloud surface. the cell center puts the cloud surface at 2/ 3 the distance between the center and the corner. This effect will occur any time a single cloud cell is present and multiple clear cells are averaged to form interpolation values. Since the interpolation results in less cloud volume, less entrainment and detrainment is measured, since VOLUME 139 there is less cloud volume with which to entrain or detrain. However, this effect only becomes an issue when a significant fraction of the cloud field has a spatial scale on the order of the model grid spacing. This can be seen by comparing vertical profiles of the total cloud core fraction for the three model resolutions (Fig. 10). The 25-m grid spacing model tetrahedral interpolation volume shows very little underestimate compared to the no-interpolation cloud core fraction, but the 100-m modelinterpolated volume is reduced in size by 50% of the uninterpolated value. This underestimate of the volume of isolated cloud cells is an inherent feature of surface interpolation, since as a cloud evaporates, the interpolated volume should decrease smoothly to zero. Thus, clouds with a volume on the order of the volume of a single model grid cell will always have their volume underestimated by interpolation. We have not been able to find a simple way to correct for this, as the effect is nonlinear and is not corrected simply by scaling the entrainment values by the ratio of the interpolated to noninterpolated cloud volumes. Therefore, the tetrahedral interpolation method must rely on having sufficient grid resolution to minimize this effect. 5. Discussion and conclusions The good agreement between our technique and the Romps method at fine grid spacing gives us confidence that our calculation is a valid representation of the flux through the cloud surface. However, as our tetrahedral interpolation redefines the cloud volume, there are interpretation issues when comparing the fluxes calculated by the two methods. We believe effects like this explain why the tetrahedral interpolation method estimates E FIG. 10. Profiles of cloud core fraction calculated without interpolation (gray line) and with tetrahedral interpolation (black line) at 25-, 50-, and 100-m grid spacings. FEBRUARY 2011 455 DAWE AND AUSTIN and D values that are smaller than those estimated via the Romps method at coarse grid resolution. Calculations of E and D with the Romps method using first-order upwind advection results in significantly larger values than when using second-order MPDATA advection. We take this to be the result of larger numerical diffusion in the upwind advection scheme. This difference is quite large, with the upwind advection giving values on the order of 50% higher than the MPDATA advection. Thus suggests that the E and D values calculated by the Romps method are very sensitive to the numerics used. Romps (2010) found that values of direct entrainment and detrainment were significantly higher than values calculated via bulk tracer budgets: our method supports this conclusion (at least at high enough grid resolution). Romps interpreted the smaller tracer budget values to be the result of assuming all entrainment and detrainment occurs with fluid having the mean properties of the environment or the cloud core, respectively. This difference between the bulk tracer and direct entrainment calculations suggests that significant recirculation occurs around the clouds, with an average air parcel entering and leaving the cloud more than once. Furthermore, Brown (1999) found entrainment and detrainment calculated via bulk tracer budgets depended on LES model resolution surprisingly little, in contrast to the strong resolution dependence shown by the direct fluxes. This suggests that changes in this recycling of cloudy air with model resolution compensate the mass flux changes resulting in little apparent E and D variability in the tracer budgets. Both these effects suggest an understanding of the dynamics of the moistened environment immediately outside the clouds, which may be important for accurate parameterization of cloud property exchanges. We have shown that by interpolating the location of cloud surfaces in an LES model we have been able to calculate reasonable entrainment and detrainment rates directly from model mass fluxes. These fluxes agree well with the direct entrainment method of Romps, with the caveat that the clouds must be well resolved by the model grid spacing so as to minimize underestimate of the cloud volume by the interpolation scheme. This tetrahedral interpolation scheme gives significant benefits over the Romps method over short time and small spatial scales, generating statistically well-behaved results suitable for use in examining entrainment and detrainment variability over the life cycle of an individual cloud. Finally, we reiterate that nothing in our technique is dependent on the shallow cumulus cloud regime, or indeed on clouds at all; it should be equally useful in calculating entrainment and detrainment through any material surface in a numerical model, subject to the caveats we have mentioned. Acknowledgments. This paper benefited from helpful discussions with Fleur Couvreux and Catherine Rio. We also gratefully acknowledge the contributions of David Romps and a second, anonymous peer reviewer, whose comments significantly improved the clarity of this paper. We thank Marat Khairoutdinov for making SAM available for our research. Support was provided by the Canadian Foundation for Climate and Atmospheric Science through the Cloud Aerosol Feedback and Climate network. All figures were generated using the matplotlib library in the Python programming language. REFERENCES Andrejczuk, M., J. M. Reisner, B. Henson, M. K. Dubey, and C. A. Jeffery, 2008: The potential impacts of pollution on a nondrizzling stratus deck: Does aerosol number matter more than type? J. Geophys. Res., 113, D19204, doi:10.1029/ 2007JD009445. Barahona, D., and A. Nenes, 2007: Parameterization of cloud droplet formation in largescale models: Including effects of entrainment. J. Geophys. Res., 112, D16206, doi:10.1029/ 2007JD008473. Bony, S., and J. Dufresne, 2005: Marine boundary layer clouds at the heart of tropical cloud feedback uncertainties in climate models. Geophys. Res. Lett., 32, L20806, doi:10.1029/ 2005GL023851. Brown, A. R., 1999: The sensitivity of large-eddy simulations of shallow cumulus convection to resolution and subgrid model. Quart. J. Roy. Meteor. Soc., 125, 469–482. ——, and Coauthors, 2002: Large-eddy simulation of the diurnal cycle of shallow cumulus convection over land. Quart. J. Roy. Meteor. Soc., 128, 1075–1093. Colman, R., 2003: A comparison of climate feedbacks in general circulation models. Climate Dyn., 20 (7–8), 865–873, doi:10.1007/ s00382-003-0310-z. Harris, J. W., and H. Stocker, 1998: Handbook of Mathematics and Computational Science. Springer-Verlag, 1056 pp. Heus, T., H. J. J. Jonker, H. E. A. V. den Akker, E. J. Griffith, M. Koutek, and F. H. Post, 2009: A statistical approach to the life cycle analysis of cumulus clouds selected in a virtual reality environment. J. Geophys. Res., 114, D06208, doi:10.1029/ 2008JD010917. Holland, J. Z., and E. M. Rasmusson, 1973: Measurement of the atmospheric mass, energy, and momentum budgets over a 500-kilometer square of tropical ocean. Mon. Wea. Rev., 101, 44–57. Khairoutdinov, M. F., and D. A. Randall, 2003: Cloud resolving modeling of the arm summer 1997 IOP: Model formulation, results, uncertainties, and sensitivities. J. Atmos. Sci., 60, 607–625. Lorensen, W. E., and H. E. Cline, 1987: Marching cubes: A high resolution 3D surface construction algorithm. Comput. Graph., 21, 163–169, doi:10.1145/37402.37422. Romps, D. M., 2010: A direct measure of entrainment. J. Atmos. Sci., 67, 1908–1927. 456 MONTHLY WEATHER REVIEW Siebesma, A. P., 1998: Shallow cumulus convection. Buoyant Convection in Geophysical Flows, E. J. Plate, Ed., Kluwer Academic Publishers, 441–486. ——, and J. W. M. Cuijpers, 1995: Evaluation of parametric assumptions for shallow cumulus convection. J. Atmos. Sci., 52, 650–666. ——, and Coauthors, 2003: A large eddy simulation intercomparison study of shallow cumulus convection. J. Atmos. Sci., 60, 1201–1219. Smolarkiewicz, P. K., and W. W. Grabowski, 1990: The multidimensional positive definite advection transport algorithm: Non-oscillatory option. J. Comput. Phys., 86, 355–375. VOLUME 139 Webb, M., and Coauthors, 2006: On the contribution of local feedback mechanisms to the range of climate sensitivity in two GCM ensembles. Climate Dyn., 27, 17–38, doi:10.1007/s00382-006-0111-2. Williams, K. D., and M. J. Webb, 2009: A quantitative performance assessment of cloud regimes in climate models. Climate Dyn., 33, 141–157, doi:10.1007/s00382-008-0443-1. Zhao, M., and P. H. Austin, 2005a: Life cycle of numerically simulated shallow cumulus clouds. Part I: Transport. J. Atmos. Sci., 62, 1269–1290. ——, and ——, 2005b: Life cycle of numerically simulated shallow cumulus clouds. Part II: Mixing dynamics. J. Atmos. Sci., 62, 1291–1310.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Faculty Research and Publications /
- Interpolation of LES cloud surfaces for use in direct...
Open Collections
UBC Faculty Research and Publications
Interpolation of LES cloud surfaces for use in direct calculations of entrainment and detrainment Dawe, Jordan T.; Austin, Philip H. Feb 1, 2011
pdf
Page Metadata
Item Metadata
Title | Interpolation of LES cloud surfaces for use in direct calculations of entrainment and detrainment |
Creator |
Dawe, Jordan T. Austin, Philip H. |
Publisher | American Meteorological Society |
Date Issued | 2011-02-01 |
Description | Direct calculations of the entrainment and detrainment of air into and out of clouds require knowledge of the relative velocity difference between the air and the cloud surface. However, a discrete numerical model grid forces the distance moved by a cloud surface over a time step to be either zero or the width of a model grid cell. Here a method for the subgrid interpolation of a cloud surface on a discrete numerical model grid is presented. This method is used to calculate entrainment and detrainment rates for a large-eddy simulation (LES) model, which are compared with rates calculated via the direct flux method of Romps. The comparison shows good agreement between the two methods as long as the model clouds are well resolved by the model grid spacing. This limitation of this technique is offset by the ability to resolve fluxes on much finer temporal and spatial scales, making it suitable for calculating entrainment and detrainment profiles for individual clouds. Copyright 2011 American Meteorological Society (AMS). Permission to use figures, tables, and brief excerpts from this work in scientific and educational works is hereby granted provided that the source is acknowledged. Any use of material in this work that is determined to be “fair use” under Section 107 of the U.S. Copyright Act or that satisfies the conditions specified in Section 108 of the U.S. Copyright Act (17 USC §108, as revised by P.L. 94-553) does not require the AMS’s permission. Republication, systematic reproduction, posting in electronic form, such as on a web site or in a searchable database, or other uses of this material, except as exempted by the above statement, requires written permission or a license from the AMS. Additional details are provided in the AMS Copyright Policy, available on the AMS Web site located at (http://www.ametsoc.org/) or from the AMS at 617-227-2425 or copyright@ametsoc.org. |
Genre |
Article |
Type |
Text |
Language | eng |
Date Available | 2011-08-01 |
Provider | Vancouver : University of British Columbia Library |
DOI | 10.14288/1.0041803 |
URI | http://hdl.handle.net/2429/32934 |
Affiliation |
Science, Faculty of Earth and Ocean Sciences, Department of |
Citation | Dawe, Jordan T., Austin, Philip H. 2011. Interpolation of LES Cloud Surfaces for Use in Direct Calculations of Entrainment and Detrainment. Monthly Weather Review, 139(2) 444–456. dx.doi.org/10.1175/2010MWR3473.1 |
Peer Review Status | Reviewed |
Scholarly Level | Faculty |
Copyright Holder | Austin, Philip H. |
Aggregated Source Repository | DSpace |
Download
- Media
- 52383-Austin_AMS_2010_2010MWR3473.pdf [ 1005.75kB ]
- Metadata
- JSON: 52383-1.0041803.json
- JSON-LD: 52383-1.0041803-ld.json
- RDF/XML (Pretty): 52383-1.0041803-rdf.xml
- RDF/JSON: 52383-1.0041803-rdf.json
- Turtle: 52383-1.0041803-turtle.txt
- N-Triples: 52383-1.0041803-rdf-ntriples.txt
- Original Record: 52383-1.0041803-source.json
- Full Text
- 52383-1.0041803-fulltext.txt
- Citation
- 52383-1.0041803.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.52383.1-0041803/manifest