<http://www.w3.org/2009/08/skos-reference/skos.html#note>
"Interpolation of LES Cloud Surfaces for Use in Direct Calculationsof Entrainment and DetrainmentJORDAN T. DAWE AND PHILIP H. AUSTINDepartment of Earth and Ocean Sciences, University of British Columbia, Vancouver, British Columbia, Canada(Manuscript received 22 April 2010, in final form 4 August 2010)ABSTRACTDirect calculations of the entrainment and detrainment of air into and out of clouds require knowledge ofthe relative velocity difference between the air and the cloud surface. However, a discrete numerical modelgridforcesthedistancemovedbyacloudsurfaceoveratimesteptobeeitherzeroorthewidthofamodelgridcell. Here a method for the subgrid interpolation of a cloud surface on a discrete numerical model grid ispresented. This method is used to calculate entrainment and detrainment rates for a large-eddy simulation(LES)model,whicharecomparedwithratescalculatedviathedirectfluxmethod ofRomps.Thecomparisonshows good agreement between the two methods as long as the model clouds are well resolved by the modelgrid spacing. This limitation of this technique is offset by the ability to resolve fluxes on much finer temporalandspatialscales,makingitsuitableforcalculatingentrainmentanddetrainmentprofilesforindividualclouds.1. IntroductionThelargestuncertaintiesinclimatesensitivityestimatesfromgeneral circulation model(GCM) simulationscomefrom the subgrid-scale parameterization of low clouds(Colman 2003; Bony and Dufresne 2005; Webb et al.2006). Specifically, stratocumulus clouds and the transi-tion regime from stratocumulus to trade cumulus are thedominant source of variance between models in the esti-mation of the cloud radiative response to changing cli-mate (Williams and Webb 2009).Propersimulationofthesubgrid-scaleeffectofcumuluscloudsinGCMsrequiresunderstandingtheratesatwhichair is entrained into and detrained from the clouds. Cloudentrainment and detrainment rates exert influences onprofiles of cloud properties, the height of the cloud tops,the amount of heat and moisture the clouds transportupward, and theheightsatwhichtheclouds depositthatheatandmoisture.Theyalsohaveeffectsontheverticaltransport of aerosols out of the boundary layer and therate at which chemical reactions can occur in thoseaerosols (Barahona and Nenes 2007; Andrejczuk et al.2008).Large-eddy simulation (LES) is an important tool usedin the study of cloud entrainment and detrainment. LESmodels apply grid resolutions on the order of 10\u2013100 m,well within the domain of the Kolmagorov 25/3 turbu-lence spectrum. This allows for a relatively simple turbu-lence model that captures the important statistics of thesubgrid-scale eddy fluxes and thus, an accurate represen-tation of the atmospheric physics in a domain \u201910 km2,whichislargeenoughtosimulateafieldofshallowclouds.LES can be assessed against results taken from field sur-veys such as the Barbados Oceanographic and Meteoro-logical Experiment (BOMEX; Holland and Rasmusson1973) or the Atmospheric Radiation Measurement Pro-gram (ARM; Brown et al. 2002), and such comparisonsshow good agreement between LES and data.Several recent studies have looked at the life cycle ofindividual clouds taken from LESs, trying to break thecloud field into its component parts (Zhao and Austin2005a,b;Heusetal.2009).Estimatesofentrainmentanddetrainment rates for individual clouds would be quiteuseful in these studies, but are difficult to achieve. En-trainment and detrainment rates for nonprecipitatingclouds are typically calculated in LES by recording bud-gets of bulk conserved tracer variables, such as the totalhumidity or the liquid water moist static energy, and in-ferring the amount of fluid exchange between cloud andclear air that is needed to explain the rate at which thattracer is being vertically advected within the cloud fieldCorresponding author address: Jordan T. Dawe, Department ofEarth and Ocean Sciences, University of British Columbia, 6339Stores Rd., Vancouver, BC V6T 1Z4, Canada.E-mail: jdawe@eos.ubc.ca444 MONTHLY WEATHER REVIEW VOLUME 139DOI: 10.1175/2010MWR3473.1C211 2011 American Meteorological Society(Siebesma and Cuijpers 1995). These budgets typically as-sumethecloudsandthecloudenvironmentarehorizontallyhomogeneousslabs;thisisamuchlessaccurateassumptionon the level of an individual cloud, and this variabilitymakesbulktracercalculationsonthisscaleproblematic.Alternatively, entrainment and detrainment could sim-plybe calculateddirectlyfromtheLESvelocity,pressure,humidity,andtemperaturefields.Siebesma(1998)definesentrainment and detrainment asE5C01A\u00FE^nC1(uC0ui),0r^n C1 (u C0 ui)dl, (1)D51A\u00FE^nC1(uC0ui).0r^n C1 (u C0 ui)dl, (2)whereEandDaretheentrainmentanddetrainmentrates(inkgm23s21), r isthedensityofair(inkgm23s21),uisthe velocity of the air (in m s21), uiis the velocity of thecloudsurface(inm s21),Aistheareaofthecloud(inm),^n is a unit vector directed out the cloud surface, and thepath integral is taken around the cloud surface at a con-stant vertical level. However, the accuracy of this methodapplied to LES suffers from the need to calculate thevelocity of the air relative to the cloud surface. In realitythese velocities are very nearly identical, ;1\u20132 m s21,butthediscretenatureoftheLESmodelgridforcesthemod-eled surface velocity to be either 0 m s21or Dx/Dt \u2019 15\u201330 m s21,whereDx is the model grid spacing andDt is themodeltimestep.Thesurfaceofthecloudonlymoveswhena grid cell\u2019s 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 av-eraging the cloud surface fluxes over the time needed foran entire grid cell to entrain or detrain. This is done bysumming the \u2018\u2018activity source\u2019\u2019 at the cloud boundary,\u203A\u203At(rA)1$C1 (ruA), (3)over the period in which the cloud surface is adjacent toagrid cell.HereA isthe \u2018\u2018activity\u2019\u2019oftheair,a scalarfieldthat is 1 inside the cloud and 0 outside. The activity canbe easily redefined to determine entrainment and de-trainment for any material surface in a numerical model;Romps usedadefinitionsimilartoSiebesmaandCuijpers(1995)\u2019s cloud core criteria (fluid with condensed liquidwater, positive buoyancy, and upward velocity) in his cal-culations. Once a grid cell is no longer adjacent to theboundary between active and inactive air, the sign of thesummed activity source is evaluated. If the summed ac-tivity source is positive, more mass has been entering theregionofactiveairthanhasbeenleaving,andthesummedactivity source is considered to be entrainment. Con-versely, a negative sum is considered to be detrainment.(In this paper, we shall refer to this method of summingthe activity source to determine entrainment and detrain-ment as the \u2018\u2018Romps method.\u2019\u2019) Romps found entrain-ment and detrainment values resulting from this methodwere approximately twice as large as those from previousbulk tracer calculations, suggesting bulk tracer calcula-tions contain significant biases. However, fine temporaland spatial resolution of entrainment and detrainmentvariability is difficult to achieve with this method becauseof the long temporal averages needed.Here we present a method for calculation of the cloudentrainment and detrainment rates that relies on inter-polation of the subgrid location of the cloud surface. Thismethod can be used to produce accurate estimates ofthe cloud entrainment and detrainment rates for indivi-dual LES clouds. In section 2 we describe this method, insection 3 we describe the model configuration we used totest this method, in section 4 we compare this calculationwith entrainment and detrainment rates calculated usingRomps (2010) direct flux calculation, and in section 5 wediscuss our results and present our conclusions.2. MethodHere we derive a method for the calculation of en-trainmentanddetrainmentfluxesthrough thesurfaceofan arbitrary volume in a numerical model. In our workthis is used to calculate the mass exchange betweencloud and environment, but this method could be usedto calculate mass flux into and out of any volume of in-terest. To maintain generality we donot yet specify howthe location of this cloud surface is determined.Consideranumericalmodelgridcellcontainingacloudsurface with normal vector C,whereC points outwardfrom the surface and has units of meters squared (Fig. 1).Thissurface,combinedwithW,theportionofthegridcellwalls that lie within the cloud (in m2), encloses a cloudvolumeVthathasunitsofmeterscubed.Siebesma(1998)givesthenetentrainmentanddetrainmentoverthecloudsurface to beE C0 D5\u00F0Cr(uiC0 u) C1 dC, (4)whereristheairdensity(inkgm23),uisthevelocityoftheair (in m s21), and uiis the velocity of the cloud interface(in m s21). Calculating this integral requires knowledge ofthe velocity field over the surface of C and the timeevolution of C, neither of which is easily calculated in anumerical model. Instead, we seek a simplified but equiv-alent calculation.FEBRUARY 2011 DAWE AND AUSTIN 445To calculate the velocity of the cloud surface, wemake use of the Leibnitz theorem:ddt\u00F0V(t)rdV 5\u00F0V(t)\u203Ar\u203AtdV1\u00F0C(t)ruiC1 dC1\u00F0W(t)ruiC1 dW. (5)Since the walls of the grid cell do not move, uiis 0 overW. If we also assume \u203Ar/\u203At \u2019 0, we getrddt\u00F0V(t)dV 5\u00F0C(t)ruiC1 dC. (6)We then combine (4) and (6) to giveE C0 D5rddt\u00F0V(t)dV C0\u00F0Cru C1 dC. (7)Next we apply the divergence theorem to simplify theflux integral through C:\u00F0V$C1 (ru)dV 5\u00F0Cru C1 dC1\u00F0Wru C1 dW. (8)As a result of the mass conservation,$C1 (ru)50, whichimplies\u00F0Cru C1 dC5C0\u00F0Wru C1 dW. (9)Substituting this into (7) results inE C0D5rdVdt1\u00F0Wru C1 dW. (10)Following Romps (2010), if (10) is positive, the resultof this calculation is assumed to all be E, and if nega-tive, D isE5max 0, rdVdt1\u00F0Wru C1 dWC18C19, (11)D5max 0, C0rdVdtC0\u00F0Wru C1 dWC18C19. (12)Therefore,wecanfindtheentrainmentanddetrainmentby calculating the rate of change of the cloud volumeinside the grid cell and the mass flux through the cloudyportion of the gridcell walls.We assume the mass flux is constant over the gridcellwall; this makes calculating mass fluxes trivial for amodelusinganArakawaC-grid,whileinterpolationcangive fluxes for other grid configurations. Expanding theterms in (10) for an Arakawa C-grid-type model with ahorizontally uniform density givesE C0 D5r(k)dV(i,j,k)dt1r(k)u(i11,j,k)Wx(i11,j,k)C0 r(k)u(i,j,k)Wx(i,j,k)1 r(k)y(i,j11,k)Wy(i,j11,k)C0 r(k)y(i,j,k)Wy(i,j,k)11/2[r(k11)1r(k)]w(i,j,k11)Wz(i,j,k11)C0 1/2[r(k)1 r(kC01)]w(i,j,k)Wz(i,j,k). (13)Thereareamultitudeofinterpolationschemesthatcanbe used to determinethe cloud volumeand cloudygridwall area in a numerical model. The simplest would beto assume that cloudy grid cells are completely filledwith cloud, and otherwise grid cells contain no cloudat all. We refer to this as the \u2018\u2018no interpolation\u2019\u2019 method.Since the cloud surface only moves when a whole gridcell undergoes condensation or evaporation this willresult in poor estimates of dV/dt, C, and W, which willresult in an overestimate of the entrainment and de-trainment. To reduce this overestimate, we interpolatethe location of the cloud surface at each time step. Nextwe describe the interpolation scheme we use for thispurpose.FIG. 1. Gridcell schematic showing geometry for the derivationof direct entrainment and detrainment fluxes using a subgrid-scaleinterpolation of the cloud surface location. Vector C is normal tothe arbitrary cloud surface, shown by the thick line. Vector W isnormaltotheportionofthegridcellwallsoccupiedbycloud,shownby the dotted lines; Wxis used to denote the x direction grid cellsurfaces; and Wydenotes the y direction surfaces. These surfacesencloseavolumeV,shownbytheshadedarea.Theuandyvectorsrepresent the locations of the horizontal velocities in the x and ydirections, respectively, for an Arakawa C-grid.446 MONTHLY WEATHER REVIEW VOLUME 139Cloud surface interpolationSeveral standard techniques exist for subgrid isosur-face interpolation in the field of computer visualization,such as the Marching Cubes algorithm (Lorensen andCline 1987). In the spirit of these techniques, we haveimplemented a scheme that relies on subdividing thegridcellsintoregularsubcells,findingthelocationofthecloud surface within these subcells, and then reassem-bling the cells to construct the location of the cloudsurface. Our interpolation scheme, which we refer to as\u2018\u2018tetrahedralinterpolation,\u2019\u2019splitsthecubicgridcellintosix pyramids, then splits each pyramid into eight tetra-hedrons (Fig. 2, left panel). This results in 48 tetrahe-drons,eachcomposedof4verticeslocatedatthegridcellcenter, the center of a gridcell wall, the center of agridcell edge, and a gridcell corner.Generally, model scalar quantities will be defined ononly one of these vertices, and so we must interpolatethe surrounding grid to determine values at the otherthreepoints.OnanArakawaC-grid,humidityisdefinedatthegridcenters,andsowemustinterpolatetofindthevaluesatthreeofthe tetrahedronvertices. Tofindvaluesat the center of the gridcell wall, we average the twovaluesonoppositesidesofthewall;atthegridcelledge,we average the four values around the edge; and at thecorners, we average the values of the eight cells sur-rounding each corner.The cloud surface itself is defined in a similar way tothe Marching Cubes algorithm (Lorensen and Cline1987). Each tetrahedron has four vertices, which caneither be inside or outside of the cloud. This results in16 possible combinations of vertex values that must beconsidered. Many of these cases share symmetries, re-ducingthenumberofindependentcasestothreeclasses.If only one of the tetrahedron vertices is cloudy (or con-versely, only one is clear) than that corner of the tetra-hedron iscut by a triangle representing the cloud surface(Fig. 2, center panel, upper left); this accounts for eightof the cases. If two vertices are cloudy and two verticesare clear, this results in the tetrahedron being cut by twotrianglesthatshareacommonedge(Fig.2,middlepanel,top right); thisaccounts for six of the cases. Finally,ifthevertices of a tetrahedron are all cloudy or clear, the sur-facedoesnotpassthroughthetetrahedron(Fig.2,middlepanel,bottom);thisaccountsfortheremainingtwocases.The locations of the vertices of the triangles that cut thetetrahedron are found by linear interpolation betweencloudy and clear tetrahedron vertices. Repeated appli-cation of this algorithm to each grid cell results in a con-tinuous surface that approximates the subgrid locationof the cloud surface (Fig. 2, right panel).Oncethegeometryofthecase isdetermined,theareaof the cloudy portion of the model gridcell wall can becalculated by dividing the cloudy area into triangles andsumming the triangle areas using the well known cross-product area formula:A5(b C0 c)3(a C0 c)jj2, (14)wherea,b,andcarethepositionvectorsoftheverticesofthe triangles. Similarly, the volume of the cut tetrahedronis calculated by subdividing it into smaller tetrahedrons,andsummingthevolumesofthesubtetrahedronsusingthetriple-product volume formula (Harris and Stocker 1998):FIG. 2. Schematic representation of the tetrahedral interpolationmethod. (left) The subdivision of the model gridcell into 48 tetrahedrons. (middle) The three cases for interpolation of the cloud surface within a subtetrahedronbetween points at the four vertices of the tetrahedron; all other tetrahedron interpolation cases are variations onthese three. Black dots represent cloudy tetrahedron points, white dots represent clear tetrahedron points, and graydotsrepresent the interpolatedcloud position betweencloudy\u2013clear tetrahedronpoint pairs.(right) A horizontal slicethrough the model\u2019s Arakawa C-grid. Rightward-pointing triangles represent u-velocity points, upward-pointingtriangles represent y-velocity points, circles represent bulk tracer quantities such as temperature or humidity, andx represents the location of the cloud surface interpolation points. Dotted lines show the boundaries of the tet-rahedral subdivision of the grid, and shading represents cloud volume.FEBRUARY 2011 DAWE AND AUSTIN 447V5(a C0 d) C1 [(b C0 d)3(c C0 d)]jj6, (15)wherea,b,c,anddarethepositionvectorsoftheverticesof the subtetrahedrons. This results in all the informa-tion needed to measure entrainment and detrainmentvia (10).3. Model descriptionWe have implemented both our tetrahedral surfaceinterpolation scheme and the time averaging scheme ofRomps (2010) in the System for Atmospheric Modeling(SAM; Khairoutdinov and Randall 2003), allowing themodel to calculate entrainment and detrainment viathese two methods simultaneously. In this paper we per-formthesecalculationsforthecloudcore,whichwedefinefollowing Siebesma and Cuijpers (1995) as fluid havingcondensed liquid water, upward velocity, and positivebuoyancy relative to the horizontal mean. The cloudcore surface then is located where the total humidityequals the saturated humidity, Tyis equal to the hori-zontal mean of Ty, and the vertical velocity is zero.Inthetetrahedralmethod,findingthelocationwhereqtequals qsat(T, p) is complicated by the nonlinearityofqsat(T, p).Forsimplicity,wedefineavariableqdiff5qt2qsat(T, p),whereqtisthetotalspecificwatercontentand qsat(T, p) is the saturated specific humidity, all inFIG. 3. Comparison of mean model profiles of core (a) entrainment and (b) detrainment calculated using thetetrahedral interpolation method (black line), the Romps method with upwind advection (dotted line), and theRomps method with MPDATA advection (gray line). (d),(e) As in (a),(b), but for the fractional entrainment anddetrainment rates.448 MONTHLY WEATHER REVIEW VOLUME 139kilograms of water per kilograms of moist air. Here qdiffcan be thought of as the moisture excess or deficiencyrelative to the saturated humidity and qdiff5 0 is the lo-cation where the total humidity equals the saturated hu-midity.Similarly,wecalculateDTy5 TyC0 Ty,wherethebar represents the horizontal mean, and Tyis the virtualtemperature in kelvin, to represent buoyancy excess ordeficiency.Verticalvelocitywrequiresnonewvariabletobecalculated,butasSAMisanArakawaC-gridmodel,itis defined at the center of the top and bottom of the gridcellinsteadofatitscenter.Therefore,wefirstinterpolatew to the center of the grid cell for consistency with lo-cation of the specific humidity and virtual temperature.Tofindthelocationofthecloudcorebetweenapointthatis in the core and a point that is outside the core, the lo-cations between these points where qdiff, DTy,andw are0 are calculated, and the value that is closest to the cloudcore point is selected as the start of the region where allthe cloud core criteria are satisfied.Since the Romps scheme relies on calculating the ad-vection of active (i.e., cloud core) air, there are actuallyseveral ways this calculation can be performed. Romps(2010)usedasimplefirst-orderupwindadvectionscheme;we have implemented both this scheme and one basedupon SAM\u2019s second-order Multidimensional PositiveDefiniteAdvectionTransportAlgorithm(MPDATA)routine (Smolarkiewicz and Grabowski 1990).We have run these schemes in a standard Global En-ergy and Water Cycle Experiment Cloud System StudyBOMEXLES(HollandandRasmusson1973;Siebesmaet al. 2003). The model was run with a domain extent of6.4 km 3 6.4 km in the horizontal and 3.2 km in thevertical. Model cloud area, vertical velocity, and cloudcore entrainment diagnosed from bulk conserved tracerbudgetsagreewellwiththeresultspresentedinSiebesmaet al. (2003) (not shown).Toexaminetheresolutiondependenceofourscheme,we ran three models: one at 25-m grid spacing in all di-rections with a 1.5-s time step, one at 50-m grid spacingwith a 3-s time step, and one at 100-m grid spacing witha 6-s time step. For most of the results we present herewe rely on the 25-m resolution model. The models wereeach run for 6 h, and the first 3 h of the simulation werediscarded as the model was still approaching the steadystate. The 15-min averages were output for the terms ofeach of our calculations. Running the 25-m resolutionmodelwiththetetrahedralsurfaceinterpolationschemeresulted in a 13% model slow down.4. Comparison with RompsHere we compare our tetrahedral interpolation methodwiththeRompsmethodcalculatedusingupwindadvectionand MPDATA advection. Comparison of the averagevalues produced by the tetrahedral interpolation methodand the Romps method with MPDATA advection over3 h of model time (Fig. 3) shows good agreement in theverticalprofilesofbothEandD;however,themagnitudeof the tetrahedral entrainment and detrainment valuesare about 15% higher than the Romps MPDATA value.The Romps method with upwind advection, however,gives results that are about 50% larger than the tetrahe-dral method. This difference is likely the result of largernumerical errors in the first-order upwind scheme com-pared with the MPDATA scheme, and is therefore anestimate of the influence of the numerical method usedon the values of E and D. In either case, the tetrahedralschemeisvastlysuperiortousingnointerpolation,whichresults in entrainment and detrainment values four timeslarger than the Romps MPDATA value (not shown).TherearesimilarlevelsofagreementbetweenC155E/Mand d 5D/M, the fraction entrainment and detrainmentrates, where M is the total vertical mass flux within thecloud core (in kg m22s21)andC15 and d have units of permeters. Since the tetrahedral method interpolation re-definestheextentofthecloudcore,wecalculateseparatevalues of M for the Romps method and the tetrahedralmethod. This results in slightly smaller M values for thetetrahedral calculation. Both C15 and d show large devia-tions between the tetrahedral and Romps calculationsFIG. 4. Average profiles of entrainment minus detrainment cal-culated using the continuity equation (dotted line), direct fluxeswith no surface interpolation (gray line), the tetrahedral interpo-lation method (thick black line), and the Romps method usingupwind advection (thin black line).FEBRUARY 2011 DAWE AND AUSTIN 449nearthecloudbaseandtop.Thesedeviationslikelyresultfrom the small cloud fractions at these heights, whichbothmakesthestatisticalsamplinglessrobustandallowssmall differences between E and D to be magnified bysmall values of M.a. Agreement with the continuity equationWhile E and D are important for calculating changesin cloud properties such as temperature and humidity,theverticalprofileofcloudfractioniscontrolledbytheirdifference,E2D.Forcomparisonwiththesedirectfluxmethods, we also calculate E 2 D from the continuityequation for a turbulent plume:r\u203Aa\u203At1\u203AM\u203Az5E C0 D, (16)where r is the density of air (in kg m23), a is the cloudfraction, and M is the fractional vertical mass flux (inkg m22s21). To satisfy continuity, the difference be-tween the amount of fluid entrained and detrained bythe clouds at a given height must equal the verticalgradientin cloudmassflux plus therate inchangeof thecloud area. We calculate the vertical gradient of massflux in (16) by interpolating vertical velocities from theedgestothecenterofeachgridcell,thenaveragingtotalmass flux within the clouds using these interpolatedvelocities. The 15-min-average values of mass flux wereoutput and the vertical derivative taken via centered dif-ferencing. Also, the 15-min-average values of \u203Aa/\u203At arecalculated within the model.Comparison of the tetrahedral values of E2Dthatwehave calculated with the Romps and no-interpolationmethods shows reasonable agreement (Fig. 4). Of thesethree calculations, the no-interpolation method is themost accurate measurement of E2D, as wecalculateitsvaluesdirectlyfromthemodel\u2019svelocityfieldwithoutanymodificationofthecloudsurface.Thisisconfirmedbytheclose agreement between the no-interpolation methodand the continuity equation. Our direct flux calculationsalsoagreewiththecontinuityequationfairlywell(Fig.4),but diverge between cloud base and 1-km height. Thetetrahedral interpolation method is biased low bythe method\u2019s redefinition of the cloud volume, while theFIG. 5. Profiles of the correlation between the Romps method and the tetrahedral interpolation method for (a) C15and (b) d. Black lines show the correlation using upwind advection in the Romps method and gray lines show thecorrelation using MPDATAadvection. The dotted line shows the 99% confidencelevel for a significant correlation.450 MONTHLY WEATHER REVIEW VOLUME 139Romps method\u2019s averaging redefines the time at whichentrainmentanddetrainmenteventstakeplace.BoththeseschemesresultinalessnegativeE2Dvaluethantheno-interpolation scheme, with the tetrahedral interpolationhaving slightly more bias. Finally, we note that theE2D,which results from the Romps MPDATA calculation(not shown), is essentially identical to the Romps up-wind calculation.b. Correlation between Romps and tetrahedralmethodsSo far we have shown that our tetrahedral interpola-tionproducesentrainmentanddetrainmentvaluesthatagreereasonablywellwiththoseproducedbytheRompsmethod. However, this does not ensure that the samevariabilityisdisplayedbythetwocalculations.Toanalyzethe variability, we take the correlation of the 15-min-averagedC15anddvaluescalculatedviathetwomethodsoverthe 3-h model run at each height to generate correlationprofiles. We do these calculations for both Romps usingupwindadvectionandRompswithMPDATAadvection.WeuseC155E/Mandd5D/MtoanalyzethevariabilityasE and D are both strongly correlated with cloud volume;more cloud volume with the same entrainment velocityresults in larger entrainment. Dividing by the mass fluxremoves this area dependence. Separate M values arecalculatedforthetetrahedralandRompsmethodstotakeinto account the volume redefinition that occurs in thetetrahedral method. Heights at which the model does nothave clouds for the entire 3-h period are excluded fromthe calculation.Correlations between the Romps method, using bothupwind and MPDATA advection, and the tetrahedralinterpolation method for both C15 and d are significant atthe 95% confidence level over most of the cloud layer(Fig. 5). The d shows nearly perfect correlation, while C15isweakerbutstillsignificant.Valuesnearthecloudbaseand top show no correlation, which we would expectFIG. 6. Time variability of (a),(c) entrainment and (b),(d) detrainment horizontally averaged over the modeldomainataheightof1 kmcalculatedvia(a),(b)thetetrahedralinterpolationmethodand(c),(d)theRompsmethodusing upwind advection.FEBRUARY 2011 DAWE AND AUSTIN 451giventhepooragreementbetweentheaveragevaluesofC15 and d between the two methods in these regions (Figs.3c,d). We take these results to show general good agree-ment between our tetrahedral interpolation scheme andthetimeaveragingschemeofRompsoverrelativelylargetime and space averages.c. Behavior over different time and space scalesThe Romps method relies on averaging over suffi-ciently long time scales for several grid cells to completean entrainment or detrainment cycle. The tetrahedral in-terpolation scheme, on the other hand, has no such lim-itation. Forexample, the tetrahedral interpolationvalueshorizontally averaged over the whole model domain ata given height show little change in variability betweeninstantaneous values, 1-min-average values, and 5-min-averagevalues(Fig.6).TheRompsmethod,ontheotherhand, showsjumps on the order of 50% of the mean fromtime step to time step. Over 1-min averages, the Rompsmethod has variability similar to that shown by the tet-rahedralinterpolationmethodwithoutanyaveraging,andFIG.7.Verticalcrosssectionofamodelcloudshowing1-min-averaged(a),(c)entrainmentand(b),(d)detrainmentcalculatedusingthe(a),(b)tetrahedralinterpolationmethodand(c),(d)theRompsmethodusingupwindadvection.Grayscale indicates amount of entrainment or detrainment and lines show the region where condensed liquid waterexists.452 MONTHLY WEATHER REVIEW VOLUME 139doesnotsettledowntoastablevalueuntil5-minaveragesare taken.When the spatial variability of the two schemes isexamined, the Romps method shows sparse measure-mentcoverage,evenwhen1-min-averagevaluesaretaken(Fig. 7). Tetrahedral interpolation values show a relativelysmooth spatial distribution at the scale of an individualcloud, and tend to be smaller than the Romps methodvalues. This is the primary advantage of the tetrahedralinterpolation method: it allows instantaneous measure-ments of entrainment and detrainment rates at the in-dividual cloud level.d. Dependence on model resolutionBoth the Romps and tetrahedral methods display rel-atively strong dependence on model resolution, with Eand D decreasing as grid resolution becomes coarser(Fig. 8). Furthermore, while E and D calculated via thetwo direct methods show good agreement at 25-m res-olution, the tetrahedral interpolation method systemat-icallyunderestimatesthevaluescalculatedbytheRompsmethod as grid spacing becomes coarser; at 100-m reso-lution, the difference in the values exceeds 50%. SincepoorestimationofmotionofthecloudsurfacerelativetotheairshouldonaverageresultinanoverpredictionofEand D, a different error must be dominating the tetra-hedral calculation.This error is likely related to a systematic under-estimate of the cloud volume as calculated by the tetra-hedral interpolation method when estimating the volumeof a cloud that occupies a singlemodel grid cell. Considersuch an isolated cloud cell, surrounded by clear grid cellsin the 2D plane (Fig. 9). For simplicity, instead of con-sidering the full cloud core calculation, let us considersimply the condensed water. Assume the cloud cell hasa qdiffvalue of 1 g kg21, and the clear cells, 21gkg21.The Romps calculation, in agreement with the modelassumptions, will treat this entire grid cell as containingcloud.The tetrahedral calculation, on the other hand, will un-derestimate the cloud volume. At the cell edge points, theqdiffvalue will be averaged between adjacent cells, givingqdiff50 and placing the cloud surface at the gridcell wall,exactly as expected. At the cell corners, however, the foursurroundingqdiffvaluesaveragetogiveqdiff520.5 g kg21,and interpolation of this value between the corner andFIG. 8. Mean (top) entrainment and (bottom) detrainment calculated at (left) 25-, (middle) 50-, and (right) 100-m grid spacing usingthe Romps method with upwind advection (gray line) and the tetrahedral interpolation method (black line).FEBRUARY 2011 DAWE AND AUSTIN 453the cell center puts the cloud surface at2/3 the distancebetweenthecenterandthecorner.Thiseffectwilloccurany time a single cloud cell is present and multiple clearcells are averaged to form interpolation values.Since the interpolation results in less cloud volume,less entrainment and detrainment is measured, sincethere is less cloud volume with which to entrain or de-train. However, this effect only becomes an issue whena significant fraction of the cloud field has a spatial scaleon the order of the model grid spacing. This can be seenby comparing vertical profiles of the total cloud corefraction for the three model resolutions (Fig. 10). The25-m grid spacing model tetrahedral interpolation vol-ume shows very little underestimate compared to theno-interpolationcloudcorefraction,butthe100-mmodel-interpolated volume is reduced in size by 50% of the un-interpolated value.This underestimate of the volume of isolated cloudcells is an inherent feature of surface interpolation, sinceas a cloud evaporates, the interpolated volume shoulddecreasesmoothlytozero.Thus,cloudswithavolumeonthe order of the volume of a single model grid cell willalways have their volume underestimated by interpola-tion.Wehavenotbeenabletofindasimplewaytocorrectfor this, as the effect is nonlinear and is not correctedsimply by scaling the entrainment values by the ratio ofthe interpolated to noninterpolated cloud volumes. There-fore, the tetrahedral interpolation method must rely onhaving sufficient grid resolution to minimize this effect.5. Discussion and conclusionsThe good agreement between our technique and theRomps method at fine grid spacing gives us confidencethat our calculation is a valid representation of the fluxthrough the cloud surface. However, as our tetrahedralinterpolation redefines the cloud volume, there are in-terpretationissueswhencomparingthefluxescalculatedby the two methods. We believe effects like this explainwhy the tetrahedral interpolation method estimates EFIG. 9. Schematic showing bias in the tetrahedral interpolationmethod,whichleadstounderestimationofcloudvolume.Rightward-pointing triangles represent u-velocity points, upward-pointing tri-angles represent y-velocity points, and circles represent bulk tracerquantity points. Black circles indicate cloudy points with total spe-cificwater1 g kg21above the saturatedspecific humidity and whitecircles indicate clear points with total specific water 1 g kg21belowthesaturatedspecifichumidity.Thegrayareaindicatestheresultingtetrahedral interpolation cloud surface.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.454 MONTHLY WEATHER REVIEW VOLUME 139and D values that are smaller than those estimated viathe Romps method at coarse grid resolution.CalculationsofEandDwiththeRompsmethodusingfirst-orderupwindadvectionresultsinsignificantlylargervalues than when using second-order MPDATA advec-tion. We take this to be the result of larger numericaldiffusionintheupwindadvectionscheme.Thisdifferenceisquitelarge,withtheupwindadvectiongivingvaluesonthe order of 50% higher than the MPDATA advection.Thus suggests that the E and D values calculated by theRomps method are very sensitive to the numerics used.Romps(2010)foundthatvaluesofdirectentrainmentand detrainment were significantly higher than valuescalculated via bulk tracer budgets: our method supportsthis conclusion (at least at high enough grid resolution).Romps interpreted the smaller tracer budget values tobe the result of assuming all entrainment and detrain-mentoccurswithfluidhavingthemeanpropertiesoftheenvironment or the cloud core, respectively. This dif-ference between the bulk tracer and direct entrainmentcalculationssuggeststhatsignificantrecirculationoccursaround the clouds, with an average air parcel enteringand leaving the cloud more than once.Furthermore, Brown (1999) found entrainment anddetrainmentcalculatedviabulktracerbudgetsdependedon LES model resolution surprisingly little, in contrasttothestrongresolutiondependenceshown bythedirectfluxes.Thissuggeststhatchangesinthisrecyclingofcloudyair with model resolution compensate the mass fluxchanges resulting inlittleapparentE and D variability inthe tracer budgets. Both these effects suggest an under-standing of the dynamics of the moistened environmentimmediately outside the clouds, which may be impor-tant for accurate parameterization of cloud propertyexchanges.We have shown that by interpolating the location ofcloud surfaces in an LES model we have been able tocalculatereasonable entrainmentanddetrainmentratesdirectly from model mass fluxes. These fluxes agreewellwiththedirectentrainmentmethodofRomps,withthe caveat that the clouds must be well resolved by themodel grid spacing so as to minimize underestimateof the cloud volume by the interpolation scheme. Thistetrahedral interpolation scheme gives significant ben-efits over the Romps method over short time and smallspatial scales, generating statistically well-behaved re-sults suitable for use in examining entrainment and de-trainment variability over the life cycle of an individualcloud.Finally, we reiterate that nothing in our technique isdependentontheshallowcumuluscloudregime,orindeedon clouds at all; it should be equally useful in calculat-ing entrainment and detrainment through any materialsurface in a numerical model, subject to the caveats wehave mentioned.Acknowledgments. This paper benefited from helpfuldiscussionswithFleurCouvreuxandCatherineRio.Wealso gratefully acknowledge the contributions of DavidRomps and a second, anonymous peer reviewer, whosecomments significantly improved the clarity of this pa-per. We thank Marat Khairoutdinov for making SAMavailable for our research. Support was provided by theCanadian Foundation for Climate and AtmosphericScience through the Cloud Aerosol Feedback and Cli-mate network. All figures were generated using thematplotlib library in the Python programming language.REFERENCESAndrejczuk, M., J. M. Reisner, B. Henson, M. K. Dubey, andC. A. Jeffery, 2008: The potential impacts of pollution ona nondrizzling stratus deck: Does aerosol number mattermorethantype?J.Geophys.Res.,113,D19204,doi:10.1029/2007JD009445.Barahona, D., and A. Nenes, 2007: Parameterization of clouddroplet formation in largescale models: Including effects ofentrainment. J. Geophys. Res., 112, D16206, doi:10.1029/2007JD008473.Bony, S., and J. Dufresne, 2005: Marine boundary layer cloudsat the heart of tropical cloud feedback uncertainties in cli-mate models. Geophys.Res.Lett.,32, L20806, doi:10.1029/2005GL023851.Brown, A. R., 1999: The sensitivity of large-eddy simulations ofshallow cumulus convection to resolution and subgrid model.Quart. J. Roy. Meteor. Soc., 125, 469\u2013482.\u2014\u2014, and Coauthors, 2002: Large-eddy simulation of the diurnalcycle of shallow cumulus convection over land. Quart. J. Roy.Meteor. Soc., 128, 1075\u20131093.Colman, R., 2003: A comparison of climate feedbacks in generalcirculationmodels.ClimateDyn.,20(7\u20138),865\u2013873,doi:10.1007/s00382-003-0310-z.Harris,J. W.,andH. Stocker,1998:HandbookofMathematicsandComputational Science. Springer-Verlag, 1056 pp.Heus, T., H. J. J. Jonker, H. E. A. V. den Akker, E. J. Griffith,M. Koutek,andF. H.Post,2009:A statistical approach to thelife cycle analysis of cumulus clouds selected in a virtual re-ality environment. J. Geophys. Res., 114, D06208, doi:10.1029/2008JD010917.Holland, J. Z., and E. M. Rasmusson, 1973: Measurement of theatmospheric mass, energy, and momentum budgets over a500-kilometer square of tropical ocean. Mon. Wea. Rev., 101,44\u201357.Khairoutdinov, M. F., and D. A. Randall, 2003: Cloud resolvingmodeling of the arm summer 1997 IOP: Model formula-tion, results, uncertainties, and sensitivities. J. Atmos. Sci.,60, 607\u2013625.Lorensen, W. E., and H. E. Cline, 1987: Marching cubes: A highresolution3Dsurfaceconstructionalgorithm. Comput.Graph.,21, 163\u2013169, doi:10.1145/37402.37422.Romps, D. M., 2010: A direct measure of entrainment. J. Atmos.Sci., 67, 1908\u20131927.FEBRUARY 2011 DAWE AND AUSTIN 455Siebesma, A. P., 1998: Shallow cumulus convection. BuoyantConvection in Geophysical Flows, E. J. Plate, Ed., KluwerAcademic Publishers, 441\u2013486.\u2014\u2014, and J. W. M. Cuijpers, 1995: Evaluation of parametric as-sumptions for shallow cumulus convection. J. Atmos. Sci., 52,650\u2013666.\u2014\u2014, and Coauthors, 2003: A large eddy simulation intercom-parison study of shallow cumulus convection. J. Atmos. Sci.,60, 1201\u20131219.Smolarkiewicz, P. K., and W. W. Grabowski, 1990: The multi-dimensional positive definite advection transport algorithm:Non-oscillatory option. J. Comput. Phys., 86, 355\u2013375.Webb,M.,andCoauthors,2006:Onthecontributionoflocalfeedbackmechanisms to the range of climate sensitivity in two GCM en-sembles. Climate Dyn., 27, 17\u201338,doi:10.1007/s00382-006-0111-2.Williams,K.D.,andM.J.Webb,2009:Aquantitativeperformanceassessment of cloud regimes in climate models. Climate Dyn.,33, 141\u2013157, doi:10.1007/s00382-008-0443-1.Zhao, M., and P. H. Austin, 2005a: Life cycle of numerically sim-ulated shallow cumulus clouds. Part I: Transport. J. Atmos.Sci., 62, 1269\u20131290.\u2014\u2014, and \u2014\u2014, 2005b: Life cycle of numerically simulated shallowcumulus clouds. Part II: Mixing dynamics. J. Atmos. Sci., 62,1291\u20131310.456 MONTHLY WEATHER REVIEW VOLUME 139"@en
<https://open.library.ubc.ca/cIRcle/collections/facultyresearchandpublications/52383/items/1.0041803>