<http://www.w3.org/2009/08/skos-reference/skos.html#note>
"Assimilating Surface Weather Observations from Complex Terrain into aHigh-Resolution Numerical Weather Prediction ModelXINGXIU DENG AND ROLAND STULLUniversity of British Columbia, Vancouver, British Columbia, Canada(Manuscript received 12 July 2005, in final form 27 June 2006)ABSTRACTAn anisotropic surface analysis method based on the mother\u2013daughter (MD) approach has been devel-oped to spread valley station observations to grid points in circuitous steep valleys. In this paper, the MDapproach is further refined to allow spreading the mountain-top observations to grid points near neigh-boring high ridges across valleys. Starting with a 3D first guess from a high-resolution mesoscale modelforecast, surface weather observations are assimilated into the boundary layer, and pseudo-upper-air data(interpolated from the coarser-resolution analyses from major operational centers) are assimilated into thefree atmosphere. Incremental analysis updating is then used to incorporate the final analysis increments (thedifference between the final analysis and the first guess) into a high-resolution numerical weather predictionmodel. The MD approaches (including one with shoreline refinement) are compared with other objectiveanalysis methods using case examples and daily mesoscale real-time forecast runs during November andDecember 2004. This study further confirms that the MD approaches outperform the other methods, andthat the shoreline refinement achieves better analysis quality than the basic MD approach. The improve-ment of mountain-top refinement over the basic MD approach increases with the percentage of mountain-top stations, which is usually low. Higher skill in predicting near-surface potential temperature is foundwhen surface information is spread upward throughout the boundary layer instead of at only the bottommodel level. The results show improved near-surface forecasts of temperature and humidity that are directlyassimilated into the model, but poorer forecasts of near-surface winds and precipitation, which are notassimilated into the model.1. IntroductionThis paper presents a method to fully incorporateboth dense local surface weather observations fromcomplex terrain and existing coarse-resolution 3Danalyses into a cost-efficient higher-resolution objectiveanalyses for local mesoscale modeling. The goal is toreduce near-surface high-resolution numerical weatherprediction (NWP) errors in mountainous terrain.Many parts of the world have complex mountainranges with valleys that are narrow and long with kinksand twists. To tackle the analysis problems in such com-plex terrain, Deng and Stull (2005, hereafter DS05) de-veloped a mesoscale analysis method using a mother\u2013daughter (MD) approach to spread valley observationsalong circuitous valleys. The MD approach is furtherrefined here for mountain-top observations, as will bedescribed in section 2.Section 3 describes the NWP model and surfaceweather observations used here. Surface observationsare valuable for mesoscale data analysis, assimilation,and forecasting (Yee and Jackson 1988; Stauffer et al.1991; Miller and Benjamin 1992; Ruggiero et al. 1996).Ruggiero et al. (2000) demonstrated that a combinationof continuous assimilation of satellite imager data andintermittent assimilation of hourly surface observationsled to a better depiction of circulations caused by cloud-shading contrasts. Alapaty et al. (2001) found that as-similating surface data into NWP models led to signifi-cant reduction in atmospheric boundary layer (BL)modeling errors, which is beneficial to subsequent airpollution forecasts.Methodology for assimilating surface observationstogether with a coarse analysis into a high-resolutionCorresponding author address: Xingxiu Deng, NumericalWeather Prediction Division, Canadian Meteorological Centre,Environment Canada, 2121 Trans-Canada Highway, Dorval, QCH9P 1J3, Canada.E-mail: xingxiu.deng@ec.gc.caMARCH 2007 DENG AND STULL 1037DOI: 10.1175/MWR3332.1\u00A9 2007 American Meteorological SocietyMWR3332mesoscale model is detailed in section 4. The method istested with case studies and daily real-time runs in sec-tions 5 and 6, respectively. A summary is provided insection 7.2. Mountain-top refinementThe MD approach in DS05 works well for surfaceobservations in valleys, but could be problematic formountain-top observations. The MD approach is notable to spread observation information from one moun-tain top to a neighboring one, because the deep valleyin between attenuates the information spread. The por-tion of information spread is called sharing factor (SF)ranging from 0 to 1. As explained in DS05, the spread-ing of observation data from a mountain top to gridpoints near other mountain tops is desirable for shallowcold-air pooling events when the mountain tops are allprotruding into the same free atmosphere air, and fordeep BL events when the mountain tops are all in thesame BL air. Therefore, the MD approach is furtherrefined here to allow such spreading by treating moun-tain-top observations differently than valley observa-tions.A modified SF for a mountain-top observation is pro-posed as follows:SoaMTH11005 SoMTH208751 H11002H20873| ZoH11002 Za|zref2H20874bH20876H208491H20850for | ZoH11002 Za| H11021 zref2; otherwise SMToaH11005 0. The super-script MT indicates mountain-top observations; thesubscripts o and a represent an observation and ananalysis grid point (GP), respectively. Here SMTois theSF at the mountain-top observation location, which is1.0. Parameter b controls the analysis decorrelationrate. The SF at any surrounding analysis GP (SMToa) de-pends upon the elevation (Z) difference between theobservation and GP, and upon the \u201Clevel-top\u201D BLdepth parameter zref2 (which was defined in DS05). Noiteration is needed.To use this approach, one must be able to distinguisha mountain-top observation from a valley observation.To do this, we first smooth the model terrain heights bythe Barnes (1964) method. The shape factor (RH9271)oftheempirical Gaussian weights is taken as 90 km, which isthe same value as used for the horizontal correlationlength scale in the first pass of the Bratseth (1986) suc-cessive-correction scheme used in the surface dataanalysis. Second, the standard deviation (H9268Z) of the dif-ference between the model terrain height and thesmoothed model terrain height is calculated within aregion of radius RH9271around each GP. If Zsis thesmoothed model terrain height, then a zeroth-order es-timate of the height of valley BL top (ZBL) can beapproximated byZBLH11005 ZsH11001 maxH208490., zref2 H11002H9268ZH20850, H208492H20850where zref2 is the level-top BL depth parameter in (1).The valley BL top here acts like a capping inversiontop, which is used to distinguish valley locations (belowthe base of the free atmosphere) from mountain-toplocations (within the free atmosphere). The approxi-mated valley BL top follows the smoothed model to-pography, but with a positive displacement (zref2 H11002H9268Z). Over complex terrain with varying ridges and val-leys, H9268Zis larger and hence the displacement is smaller.Over flat terrain, H9268Zis smaller and therefore the dis-placement is larger. In the case of zero H9268Z, the displace-ment equals the level-top BL depth parameter zref2used in the SF equations [(1) in this paper and (5) inDS05].If the unsmoothed terrain height for any GP is aboveZBL, then that GP is assumed to be in the free atmo-sphere, and is treated as a mountain-top location. Theapproximation of ZBLis crude, but provides a simpleand effective way to distinguish mountain-top from val-ley observations. As mentioned in DS05, for any stationnot collocated with a model GP, the station is first ap-proximated as being collocated at the nearest-neighboring model GP having the least elevation dif-ference between them. In such a case, the elevation ofthe nearest-neighboring model GP is compared withZBL. Any surface weather station that is above ZBListreated as a mountain-top station, whereas any surfacestation that is below ZBLis treated as a valley station.The SF for a mountain-top station is obtained through(1), whereas the SF for a valley station is calculated viaEq. (5) in DS05.To illustrate, Fig. 1 shows the model topography oftest domains of 3- and 2-km horizontal grid spacing insouth British Columbia (BC), Canada, where our as-similation and forecast experiments are conducted.Figure 2 shows cross sections of the model topogra-phy, the highly smoothed model topography (Zs), theapproximated height of valley BL top (ZBL), and H9268Z.The top (bottom) panel in Fig. 2 is along the projectionline in the top (bottom) panel of Fig. 1. Any surfacestation or GP above the dot\u2013dashed ZBLline will betreated as a mountain-top station or GP.To further illustrate, Fig. 3 shows terrain heights inthe Coast Mountains of BC, just north of Vancouver[Lower Fraser Valley (LFV) in Fig. 1]. With the samecase example (7\u20138 March 2003) as in DS05, a virtualmountain-top observation at o3 (Fig. 3) is extractedfrom the \u201Ctruth\u201D model surface forecasts of \u201Cfraternaltwin\u201D experiments (see DS05). The observation incre-1038 MONTHLY WEATHER REVIEW VOLUME 135ment (observation minus the first guess from the\u201Canalysis\u201D model forecasts) is H110021.56 K. The analysisincrements (analysis minus the first guess) produced bythe observation o3 using the MD approach after moun-tain-top refinement (method MD_MT) are shown inFig. 4. Method MD_MT successfully spreads the moun-tain-top observation increment to other mountain tops,while giving very small spread in the valleys.In DS05, the MD approach was compared with twoother published methods: GAUSS is the generic Gauss-ian spread built into the Advanced Regional PredictionSystem (ARPS; Xue et al. 2000) Data Assimilation Sys-tem (ADAS; Brewster 1996); and TERR_DIFF is theterrain difference method of Miller and Benjamin(1992). See DS05 for more details about GAUSS andTERR_DIFF. Redoing that fraternal twin experimentfrom DS05, but now including method MD_MT, yieldsthe analysis verification statistics in Tables 1 and 2.When only the mountain-top observation at o3 is ana-lyzed, Table 1 shows that MD_MT verifies second bestFIG. 1. The model topography (darker shading indicates higherelevations) for (top) 3- and (bottom) 2-km grid spacing. The blacklines are the projection lines for the cross sections in Fig. 2. SeeFig. 5 for the geographic relationship between these two domains.The cities of Vancouver and Victoria, British Columbia, are alsoindicated with unlabeled circles, near the LFV and the southerntip of Vancouver Island, respectively.FIG. 2. Cross sections of the model topography (solid line), thesmoothed model topography (dashed line), the approximatedheight of valley BL top (dot\u2013dashed line) and the std dev of thedifference between the model topography and the smoothedmodel topography (dotted line). (top) The cross section is along49.0\u00B0N cutting through the 3-km domain (see the black line in thetop of Fig. 1). (bottom) The cross section is along 50.08\u00B0N cuttingthrough the 2-km domain (see the black line in the bottom of Fig.1). Any surface station located above the dot\u2013dashed line istreated as a mountain-top station. The same is true for GPs;namely, the regions where the solid line is above the dot\u2013dashedline are considered to be mountain-top locations.MARCH 2007 DENG AND STULL 1039(after MD) against valley stations only, and best whenverified against mountain stations only. When two val-ley observations (o1 and o2 from Fig. 3) and one moun-tain-top observation (o3) are analyzed, Table 2 showsthat MD_MT verifies best against all the other valleyand mountain stations combined.To further test MD_MT, we next describe the NWPmodel and methodology used for one case study (a for-est fire event) and for 2 months of daily real-time runsusing real data.3. The numerical model and dataa. The numerical modelThe NWP model into which the surface data are as-similated is the Mesoscale Compressible Community(MC2) model (Laprise et al. 1997; Benoit et al. 1997)version 4.9.1. It uses vertically stretched Gal\u2013Chen ter-rain-following coordinates to obtain greater resolutionnear the surface.MC2 is configured with five one-way self-nested gridswith horizontal grid spacings of 108, 36, 12, 4, and 2 or3 km. The National Centers for Environmental Predic-tion (NCEP) Eta Model1analysis and forecasts fromthe \u201C104\u201D grid (with a 90.7-km grid spacing true at60\u00B0N) were used as the initial and boundary conditionsfor the coarsest MC2 grid. The Eta forecasts were avail-able every3hupto84h.Westart from the 104-gridoutput because this domain extends far enough westover the Pacific to reduce upstream boundary errors,and also because this output was continuously availableand timely accessible via the Internet.The 4-, 2-, and 3-km meshes (Fig. 5) have 35 layers(18 below 1500 m AGL) in the vertical, with the modeltop at 23 km. The first \u201Cthermodynamic\u201D (\u201Cmomen-tum\u201D) level is located at 5.3 (10.6) m above the modelground.Objective analyses of surface potential temperatureand specific humidity are performed at the first ther-modynamic level. A data assimilation (DA) experimentfor the case study is performed on the 3-km domain(Fig. 5 and Fig. 1, top), centered at Vernon in south-1Effective 25 January 2005, the Eta Model has been renamedthe North American Mesoscale (NAM) model.FIG. 4. The AIs (isopleths) from method MD_MT for potentialtemperature at the lowest (terrain following) model level in re-sponse to a single mountain-top observation at o3 (whitesolid triangle). The contour interval is 0.2 K. Darker shading in-dicates higher elevations (m). The observation increment at o3 isH110021.56 K.FIG. 3. Virtual surface observation stations (solid triangles)used in the analysis and verification. The stations are positioned atthe truth model terrain height in the Coast Mountains north ofVancouver. Station o1 (50.326\u00B0N, 123.578\u00B0W) is near the mouth ofthe Elaho River. Station o2 (50.3344\u00B0N, 122.767\u00B0W) is near thetown of Pemberton. Station o3 (50.3773\u00B0N, 123.2363\u00B0W) is overthe mountain top. Terrain elevations (m) are from the truth model(MC2 2 km). Darker shading indicates higher elevations, with amaximum elevation difference of 2055 m in this figure. Thedashed line in lower-left corner is a fjord (Jervis Inlet).TABLE 1. The NRMSEs between the analyses and observationsfor different sets of the verification stations. One set is the valleyobservations (b1, b2, b3, and b4) from Fig. 3; the other is themountain-top observations (m1, m2, and m3). The analyses wereproduced when only the single mountain-top observation o3 isused. The observation stations are shown in Fig. 3. NRMSE isRMSE for each method normalized by the RMSE of the FSTG.Smaller NRMSE corresponds to better analyses. NRMSE close to1.0 indicates that the observations contributed relatively littlevalue relative to the first guess.MethodVerified againstb1, b2, b3, and b4Verified againstm1, m2, and m3FSTG 1.0 1.0GAUSS 1.103 0.275TERR_DIFF 1.406 0.214MD 1.001 1.0MD_MT 1.003 0.1631040 MONTHLY WEATHER REVIEW VOLUME 135central BC. Real-time DA runs are performed for the2-km mesh (Fig. 5 and Fig. 1, bottom). For all DA runson 3- or 2-km meshes, five GPs are cut from each sideof the analysis domains shown in Fig. 5, to reduce theeffects of lateral boundary errors.b. DataHourly surface weather observations are from the real-time \u201CEmergency Weather Net Canada\u201D (EmWxNet)database, which includes surface temperature, relativehumidity, wind speed/direction, and mean sea levelpressure (or surface pressure at some stations). Tem-perature and humidity are used for both assimilationand verification, while other variables are used for fore-cast verification only.Potential temperature varies smoothly over moun-tainous terrain when the BL is relatively deep and wellmixed (Miller and Benjamin 1992); thus, it is chosen asthe temperature variable for analysis. It is also the tem-perature variable analyzed in ADAS.Specific humidity is chosen as the moisture variablefor analysis and assimilation into MC2 for two reasons:1) MC2 uses specific humidity; and 2) specific humidityis a continuous variable. Specific humidity is also themoisture variable analyzed in ADAS. To derive spe-cific humidity from relative humidity, the saturationspecific humidity is first calculated using the enhancedTeten\u2019s formula, as in ARPS. This calculation requiressurface temperature and station pressure as input.There are few reports of station pressure from theEmWxNet. In the case of missing station pressure,mean sea level pressure is used to approximate stationpressure according to the altimeter equation (standardatmosphere) for stations with elevations below 500 m.4. MethodologyAssimilating surface weather observations into ahigh-resolution NWP model in this study involves threemajor components: horizontal spreading, verticalspreading, and data insertion.a. Horizontal spreading of surface observationsThe isotropic assumption typically used in back-ground error correlation models for horizontallyspreading observations is not valid for mountainous ter-rain. Instead, surface observations are analyzed aniso-tropically using MD_MT with ADAS to generate a sur-face data analysis (see DS05). This is done by optimallycombining surface observations and a high-resolutiongridded first guess from the lowest terrain-followinglevel of a previous MC2 run. The parameters used ineach analysis method are the same as in DS05.b. Vertical spreading of surface observations andcombination with pseudo-upper-air dataSurface observations are available at only one ter-rain-following level. Meanwhile, major operational cen-ters generate 3D meteorological analyses daily by as-similating many types of measurements (rawinsondes,aircraft, wind profilers, satellite and radar, etc). Thissection describes how to spread surface information up-ward, and to combine it with pseudo-upper-air datainterpolated from a coarser, 3D analysis from the majoroperational centers.FIG. 5. MC2 grid domains for H9004x H11005 4, 3, and 2 km. MC2 4-kmoutput provides nesting files to initialize 2- and 3-km runs. MC2 at3 or 2 km is run to provide first-guess fields for analysis, and alsofor assimilation runs after analysis. Five grid points are cut fromeach side of the analysis domain shown here, to reduce lateralboundary effects in the assimilation. The forest fire location nearthe town of McLure is indicated with the X.TABLE 2. Same as in Table 1, but for all of the verificationstations (the valley stations b1, b2, b3, and b4, and the mountain-top stations m1, m2, and m3 from Fig. 3). The analyses wereproduced using two observations (o1 and o2) in different valleysand one mountain-top observation o3.MethodVerified against b1, b2, b3, b4,m1, m2, and m3FSTG 1.0GAUSS 1.759TERR_DIFF 1.967MD 0.972MD_MT 0.415MARCH 2007 DENG AND STULL 1041To obtain pseudo-upper-air data on the terrain-following MC2 3-km (or 2 km) model levels, MC2 isutilized as an interpolation tool here. This is done byintegrating MC2 for only 1 h and self-nesting the coars-est grid (108 km) successively down to 36, 12, 4, andfinally to the finest one (3 or 2 km). The resulting MC23-km (or 2 km) output on terrain-following levels at 0 hare the desired pseudo-upper-air data. By assimilatingthis interpolated Eta data, we effectively incorporate allthe satellite, radar, aircraft, and rawinsonde data thatwere assimilated by NCEP.At this stage, we have a surface data analysis,pseudo-upper-air data, and the first guess from a pre-vious MC2 3-km (or 2 km) forecast, all on the MC2model levels and at every model grid column. In devis-ing a profile (PROF) scheme to combine the three setsof data, we adopt assumptions similar to the ones ofYee and Jackson (1988): surface observations describeatmospheric state in the BL, while coarser, 3D analysesfrom major operational centers provide informationabove the BL. In the PROF scheme, potential tempera-ture and specific humidity analyses at the lowest modellevel are applied to the whole BL. Above the BL top inthe free atmosphere, the final analysis FAis a weightedaverage of the first guess FBand the pseudo-upper-airdata FUA:FAH20849r, kH20850H11005 H208511 H11002 WH20849r, kH20850H20852FBH20849r, kH20850H11001 WH20849r, kH20850FUAH20849r, kH20850,H208493H20850where indexes r and k correspond to horizontal positionand vertical level, respectively. Weights W(r, k) in (3)depends on the ratio of error variances of the pseudo-upper-air data and first guess:WH20849r, kH20850H1100511 H11001H20875H9268UAH20849r, kH20850H9268BH20849r, kH20850H208762. H208494H20850As the pseudo-upper-air data are obtained by inter-polating the Eta analysis to the MC2 3-km (or 2 km)model levels, the error variance of the pseudo-upper-air data is assumed to be the same as that of the firstguess: H9268UA(r, k) H11005 H9268B(r, k).The BL height is diagnosed from a profile methodusing a slab idealization of the mixed layer (Stull 2000),where the potential temperature jump in the entrain-ment zone is taken as H9004H9258H11005 1.5 K. Namely, the analyzedpotential temperature at the lowest model level H9258SA(1)is compared with the potential temperature profile ofthe pseudo-upper-air data (H9258UA) at successively highergrid points. When H9258SA(1) H11002 H9258UA(1) H11022H110021.5 K, the BLdepth Ziis the height (above the model ground) ofthe model level k at which the criterion of H9258UA(k) H11002H9258SA(1) H11350 1.5 K is first met. When H9258SA(1) H11002 H9258UA(1) H113491.5 K, Ziis assumed to be not shallower than 300 m,which represents nocturnal BLs or shallow BLs withcold-air pooling. Figure 6 illustrates the relationship be-tween the BL depth Ziand H9258SA(1) H11002 H9258UA(1).c. Insertion of analysis incrementsThe objective analysis procedures used in this work(same as in DS05) do not provide mass and motionfields that are in optimal balance to initiate NWP. Thisstudy uses the incremental analysis updating (IAU;Bloom et al. 1996), instead of an explicit initializationprocedure. In IAU, the analysis increments (AIs) aregradually incorporated into the MC2 model. BecauseMC2 uses a process-splitting method (Bergeron et al.1994) for its different prognostic terms, it was conve-nient to apply the AIs after the physical parameteriza-tions in the vertical, but before the treatment of hori-zontal diffusion and the Asselin\u2013Robert time filter.IAU is implemented in MC2 largely following theexisting ADAS IAU (Brewster 2003; Bloom et al.1996). Parallelization of IAU uses the software archi-tecture of MC2. In their applications of ADAS IAU(Brewster 2003; Xue et al. 2002), the analysis incre-FIG. 6. Schematic diagram illustrating the relationship of the BLdepth Ziand the difference between the analyzed potential tem-perature and the interpolated Eta analysis at the bottom modellevel [H9258SA(1) H11002 H9258UA(1)]. The actual slope of the diagonal linedepends on the upper-air profile.1042 MONTHLY WEATHER REVIEW VOLUME 135ments for radar data were introduced to a storm-scalenumerical model employing a constant time weightingover a 10-min period. In our work, the analysis incre-ments are introduced into the model over a 1-h IAUwindow using a constant time weighting (sections 5and 6).5. Case study resultsNumerical experiments are performed for 29\u201330 July2003 when a large forest fire occurred near McLure, BC(see Fig. 5). During that time, fair weather associatedwith a 50-kPa ridge over southern BC allowed BL pro-cesses to dominate. The 3-km domain (Fig. 5) is usedfor this assimilation case study.The Eta Model analysis and forecasts, started at 0000UTC 29 July 2003, are used to drive the 108-km grid,which in turn drives the nests of finer grids (see Fig. 7).The MC2 3-km control run (CTRL) is started at 1200UTC 29 July 2003 from the MC2 4-km output, andprovides its 12-h output as the first guess for the analy-sis at 0000 UTC 30 July 2003. Surface analyses are per-formed for 0000 UTC 30 July 2003 by blending hourlysurface observations with the first guess at the lowestterrain-following model level valid at the same time.Those surface stations with a difference between theiractual elevation and model topography greater than500 m are excluded from analysis and verification. Outof the total 134 stations in the DA domain, 20 stationswere excluded for this reason.Pseudo-upper-air data are obtained from the Etaanalysis valid at 0000 UTC 30 July 2003. A final high-resolution analysis at 0000 UTC 30 July 2003 is pro-duced by combining the surface analysis, pseudo-upper-air data, and first guess using the PROF scheme.A DA run at 3-km grid spacing is started at 0000 UTC30 July 2003 from the final analysis via the IAU tech-nique. Lateral boundary conditions for the DA runsand for the CTRL run are the same, and are from theMC2 4-km output. Verification of subsequent modelforecasts at the lowest model level against surface ob-servations is performed during a 24-h forecast periodfor the DA forecast run, and are compared with veri-fication of the CTRL run during the 12\u201336-h forecastperiod from 0000 UTC 30 July to 0000 UTC 31 July2003. The verification of both the surface analysis andforecast fields are presented next.Note that the comparison against this CTRL run isnot strictly optimal, since this CTRL run starts earlierthan the DA runs. This suboptimal choice of CTRLimplies a probable favor in DA experiments comparedto CTRL. A much more appropriate CTRL run wouldstart from the analysis of an operational model valid atthe same initial time as that of all other experiments.The purpose of this paper is not only to demonstratethe usefulness of surface data but also to assess theimpact of data assimilation on real-time high-resolutionmodel runs that do not have DA module originally. Inthis paper, we sacrifice the optimum, in order to verifythe assimilation strategy over a long time period in realtime without modifying the current daily real-time fore-cast system. Nevertheless, one should keep this prob-lem of suboptimal CTRL in mind when comparing theresult from one DA experiment against that fromCTRL. The problem with suboptimal CTRL is mini-mized when comparing the results from two DA ex-periments.a. Analysis verificationFigure 8 shows all 118 surface stations within theanalysis domain. The domain for DA runs is smallerthan the analysis domain, leaving 114 stations withinthe DA domain for verification. Among them, 65 sta-tions (indicated by closed triangles) reported tempera-ture observations at 0000 UTC 30 July 2003. Observa-tions separated by 100 m or less in the horizontal andvertical are averaged to create a smaller number of\u201Csuperobservations.\u201D Because only a subset of thesestations also reported relative humidity and pressure,there are only 16 stations for which specific humiditycould be calculated (Fig. 9).Table 3 shows that all methods give improved poten-tial temperature analyses compared to the first guess(FSTG) as measured by bias, mean absolute error(MAE), and normalized RMSE (NRMSE). Both MDand MD_MT outperform GAUSS and TERR_DIFF;MD produces the lowest bias and MAE, whereasFIG. 7. Schematic diagram illustrating model time lines of theMC2 self-nested grids including the 3-km CTRL and DA runs.MC2 4-km provides boundary conditions for the 3-km CTRL runand the 3-km DA runs. SFC indicates surface data.MARCH 2007 DENG AND STULL 1043MD_MT has the lowest NRMSE. The improvement ofMD_MT over MD is not as large as that for the virtualobservation case in section 2. One reason is that only5% of the stations (3 out of 65) are identified as moun-tain-top observations in this case, while 33% (1 out of3) were from mountain top for the virtual observationcase.For specific humidity, Table 4 shows that MD andMD_MT outperform GAUSS and TERR_DIFF. Nospecific humidity observations are from the mountain-top stations, therefore MD_MT H11005 MD. The MD ap-proach with shoreline refinement (MD_LSMG, seeDS05) is not studied here, because there are no stationsover water in this domain.The results above confirm the conclusion from DS05:the MD approach outperforms GAUSS and TERR_DIFF in producing better surface analyses. The resultsalso show that the mountain-top refinement provides asmall improvement over the basic MD approach. Insection 5b, the impacts of assimilating surface weatherobservations using MD_MT on subsequent near-sur-face forecasts are examined.b. Verification of MC2 forecasts1) FORECAST EXPERIMENTSEight experiments are conducted: one control experi-ment without assimilation (CTRL) and seven other ex-periments, which assimilate various combinations ofthe surface and pseudo-upper-air data, use different in-sertion rates, or assimilate different meteorologicalfields (see Table 5).Experiment SURFDA assimilates the surface datainto only the lowest model level, and uses the first guessat all higher model levels. Experiment SF_PROF usesthe surface analysis only, but the AIs at the lowestmodel level are spread vertically by the PROF schemeto the whole BL. Experiment PROFDA uses the sur-face and pseudo-upper-air data combined by the PROFscheme. By comparing SF_PROF and PROFDA, onecan see if the assimilation of surface data plays a dom-inant role on subsequent near-surface forecasts. For theabove three experiments, the temperature AIs are in-corporated all at once within a single time step (30 s)window.FIG. 9. Surface station locations (closed triangles) with suffi-cient data to calculate specific humidity at the analysis time (0000UTC 30 Jul 2003) and the model topography (m). Darker shadingindicates higher elevations.TABLE 3. Verification of analyzed potential temperatures for allreporting verification stations at the analysis time: 0000 UTC 30Jul 2003. Here N equals the total number of observations. For allthese statistics, smaller is better.Method N Bias (K) MAE (K) NRMSEFSTG 64 H110025.258 7.687 1.0GAUSS 64 0.758 2.099 0.515TERR_DIFF 64 0.712 2.324 0.546MD 64 0.089 1.600 0.391MD_MT 64 0.133 1.610 0.387FIG. 8. Surface weather stations (triangles) are superposed onthe model topography (m) for the 3-km domain shown in Fig. 5.Darker shading indicates higher elevations. Open triangles repre-sent station locations with missing temperature observations atthe analysis time (0000 UTC 30 Jul 2003), whereas closed trianglesindicate station locations reporting temperatures at the analysistime.1044 MONTHLY WEATHER REVIEW VOLUME 135Experiments TH1 (i.e., experiment PROFDA)through TH6 assimilate surface and upper-air tempera-ture fields employing different insertion rates over dif-ferent IAU time windows. TH1 incorporates AIs all atonce within a single 30-s time step. Over a 1-h window,experiments TH2, TH3, and TH6 incorporate AIs ev-ery 30, 120, and 1200 s, respectively. The 1-h window isselected because the BL responds to surface forcingswith a time scale of about1horless (Stull 1988).The settings for experiments ATQ and AQ are thesame as those for experiment TH6, except that ATQassimilates both temperature and specific humidity, andAQ assimilates only specific humidity.2) VERIFICATION RESULTSStatistical assessments during the 1\u201312-h forecast pe-riod in Table 6 suggest that different DA experimentshave variable success at predicting near-surface fields.The seven DA experiments can be divided into threegroups. Comparisons are made below within each group.For potential temperature forecasts, experimentsSURFDA, SF_PROF, and PROFDA outperformCTRL. By assimilating surface temperature at only thelowest model level, SURFDA produces very little im-provement over CTRL. When surface potential tem-perature analyses are spread throughout the whole BL(viz. experiment SF_PROF), larger improvement forpredicting surface potential temperatures are achieved.By combining surface and pseudo-upper-air data,PROFDA gains only slightly larger improvement thanSF_PROF. Thus, assimilating surface data that werespread vertically over the BL was valuable for this for-est fire weather case.Verification of mean sea level pressure (SLP) fore-casts reveals similar performance. All experiments ex-cept SURFDA outperform CTRL. Overall, improve-ment of SLP forecasts is smaller than that of potentialtemperature. This is reasonable as temperature, notpressure, observations are directly assimilated into themodel. For wind forecasts, all three DA experimentsunderperform CTRL in terms of normalized root-mean-square vector error (NRMSVE) due to initial im-balances between mass and wind fields caused by asudden change in temperature fields. SURFDA has thesmallest NRMSVE, probably because the assimilatedtemperature information is lost soon after insertion.The impact of varying insertion rate on potentialtemperature forecasts is much larger than on SLP andwind forecasts (cf. experiments TH1, TH2, TH3, andTH6 in Table 6). The errors of potential temperaturedecrease significantly from experiment TH1 to TH2 (orTH3 and TH6). This implies that introducing the AIsover a 1-h window helps the model to better retain theassimilated information than adding the increments atonly the initial time. The differences between the errormeasures of experiments TH2, TH3, and TH6 aresmall. Comparatively, TH6 has the smallest errors forpotential temperature and vector wind.The top panel of Fig. 10 shows NRMSE of potentialtemperature versus forecast hour for experimentsCTRL, TH1, TH3, and TH6. The NRMSEs for the DAexperiments are much smaller than 1.0 (NRMSE ofCTRL) during the first 12 h (except at 4 h, see expla-nation for the exception later in this section), and laterincrease gradually. Finally the NRMSE becomes closeto 1.0. By applying the AIs over a 1-h window ratherthan all at once, TH3 and TH6 produce much smallerNRMSE than TH1 for the first 7 h (except 4 h) offorecasts. Reducing the insertion rate at which the AIsTABLE 5. Experimental design. Variable T is temperature and q is specific humidity.Expt name Surface analysis UA analysis Vertical spreading DA window Insertion rate Variables assimilatedCTRL NoSURFDA Yes No No 30 s 30 s TSF_PROF Yes No Yes 30 s 30 s TPROFDA (TH1) Yes Yes Yes 30 s 30 s TTH2 Yes Yes Yes 1 h 30 s TTH3 Yes Yes Yes 1 h 120 s TTH6 (AT) Yes Yes Yes 1 h 1200 s TAQ Yes Yes Yes 1 h 1200 s qATQ Yes Yes Yes 1 h 1200 s T, qTABLE 4. Same as in Table 3, but for verification of analyzedspecific humidity.Method NBias (1.0 H1100310H110024kg kgH110021)MAE (1.0 H1100310H110024kg kgH110021) NRMSEFSTG 16 14.941 21.921 1.0GAUSS 16 H110023.595 7.878 0.621TERR_DIFF 16 H110024.480 9.440 0.665MD 16 H110023.181 6.770 0.594MD_MT 16 H110023.181 6.770 0.594MARCH 2007 DENG AND STULL 1045are introduced decreases NRMSE mainly for the first 1h of forecasts. Reducing the insertion rate over a 1-hwindow implies a larger fraction of the AIs is applied ateach insertion step and thus experiences less diffusionafter 1 h.The time series of NRMSVE for experiments CTRL,TH1, TH3, and TH6 are shown in the bottom panel ofFig. 10. By assimilating temperatures only, the modelproduces poorer wind forecasts than CTRL during thefirst 12 h due to initial imbalances between the massand wind fields. During the second 12 h, all DA experi-ments and CTRL produce similar forecasts. This couldbe because the assimilated information propagates outof the domain. Another reason could be that the massand wind fields are adjusted to be in balance. By ap-plying the temperature AIs over a 1-h window ratherthan all at once (thus reducing initial imbalances), TH3and TH6 decrease NRMSVE for the 1-h forecast com-pared to TH1, but increase NRMSVE a little during the2\u20138 h. Reducing the rate at which the AIs are intro-duced further decreases NRMSVE for the 1-h forecastwhile reducing the growth of NRMSVE later.The summed total precipitation amounts over the en-tire domain for experiments CTRL, TH1, TH2, TH3,and TH6 are listed in Table 7. Actual observations (notlisted) show trace to very light precipitation only atsome stations during the simulation period. Comparedto CTRL, data assimilation experiments result in muchmore precipitation during the first 3 forecast hours, andslightly smaller precipitation afterward. By applying theAIs over a 1-h IAU window rather than all at once,TH2, TH3, and TH6 reduce excessive precipitation dur-ing the first 3 forecast hours compared to TH1. But theprecipitation forecasts from the DA runs remain poorerthan from the CTRL run. This may imply that a betterinitialization step is required when assimilating surfacedata into the NWP model. Also, in steep mountainousregions much of the precipitation is orographicallymodified, so poor wind forecasts lead to poor precipi-tation forecasts.FIG. 10. Time series, for (top) NRMSE of surface potentialtemperature and (bottom) NRMSVE of surface winds from ex-periments CTRL, TH1, TH3, and TH6. Smaller NRMSE (orNRMSVE) corresponds to better forecasts.TABLE 6. Verification of surface potential temperature (H9258), vector wind (V), and MSLP forecasts in terms of bias during a 12-hforecast period from 0100 to 1200 UTC 30 Jul 2003. NRMSE (NRMSVE) is RMSE (RMSVE) for each experiment normalized by theRMSE (RMSVE) of the control run (CTRL). The DA experiments are defined in Table 5. For all these statistics, smaller (closer tozero) is better.1\u201312-h forecastH9258(K) (N H11005 601) SLP (kPa) (N H11005 153) V (m sH110021)(N H11005 506)Expt Bias MAE RMSE NRMSE Bias MAE NRMSE RMSVE NRMSVECTRL H110024.838 5.666 6.472 1.000 0.376 0.397 1.000 1.745 1.0SURFDA H110024.853 5.623 6.442 0.996 0.376 0.397 1.000 1.750 1.003SF_PROF H110024.328 5.288 6.089 0.941 0.321 0.365 0.951 1.968 1.128PROFDA H110024.163 5.179 5.979 0.924 0.301 0.357 0.935 1.996 1.144TH2 H110023.846 4.885 5.679 0.878 0.269 0.363 0.940 2.024 1.160TH3 H110023.871 4.910 5.700 0.881 0.270 0.362 0.939 2.023 1.160TH6 (AT) H110023.752 4.792 5.600 0.865 0.277 0.358 0.936 1.991 1.141AQ H110024.921 5.731 6.554 1.013 0.385 0.408 1.021 1.773 1.016ATQ H110023.688 4.742 5.548 0.857 0.275 0.357 0.924 2.033 1.1661046 MONTHLY WEATHER REVIEW VOLUME 135Now we look at the impact of assimilating differentmeteorological fields on subsequent model forecasts.For potential temperature and SLP forecasts, experi-ment AQ (assimilating specific humidity only) pro-duces poorer forecasts than CTRL (Table 6). The fore-casts from experiments AT and ATQ are greatly im-proved compared to the CTRL forecast. Byassimilating both temperature and specific humidity,ATQ is improved compared to AT that assimilatestemperature only.The DA experiments AT, AQ, and ATQ give poorerwind forecasts than CTRL. In contrast to the errormeasures of potential temperature and SLP, NRMSVEfor wind is the smallest from AQ, and the largest fromATQ. Winds are more sensitive to temperature pertur-bations than to moisture perturbations.Verification statistics of specific humidity for experi-ments AT, AQ, and ATQ are shown in Table 8. As-similating only specific humidity (experiment AQ)greatly decreases bias, MAE, and NRMSE in the fore-casts of specific humidity compared to CTRL. ATQ hassmaller MAE and NRMSE than AT.The impact of assimilating different meteorologicalfields on the diurnal variation of surface temperature isexamined at three stations near McLure and three sta-tions in the Okanagan Valley (Fig. 11). McLure and theOkanagan Valley suffered extensive forest fires in thepast few years, for which improved forecasts could haveaided fire fighters in saving more homes. The diurnalvariations of surface temperature forecasts for experi-ments AT and AQ are almost the same as experimentsATQ and CTRL, respectively, and therefore not in-cluded here.Figure 12 compares near-surface temperature fore-casts (at the lowest model level, which is 5.3 m abovemodel ground) from experiments CTRL and ATQ tothe observed surface temperatures at three stationsnear McLure every 1 h during the 24-h forecast period.The first guess (from CTRL) underforecasts surfacetemperature at the analysis time (0000 UTC 30 July2003). By assimilating local surface temperature obser-vations, experiment ATQ produces improved forecasts.The improvement decreases from 1- to 2-h forecast forthese three stations.Time series of observed and forecasted surface tem-peratures are also shown for three stations located inthe populated Okanagan Valley in Fig. 13, where simi-lar conclusions can be drawn. However, the differencein the magnitude between the observed and forecastedsurface temperature is obvious, caused partly by thedifference between actual and model elevation (Table9) in steep terrain.At 0400 UTC, Figs. 12 and 13 show either a largedrop in surface temperature at some stations or missingreports at others. This helps explain the sudden jumpof NRMSE for surface potential temperature at4hinFig. 10.FIG. 11. Three surface stations (closed triangles) near McLure(a forest fire location indicated by a closed square) and threeothers farther south in the Okanagan Valley. The model terrainheights (m) are shown by shading. Darker shading corresponds tohigher elevations. The shaded region corresponds to the 3-kmdomain in Fig. 5.TABLE 8. Verification of surface specific humidity (q) for the1\u201312-h forecast period from 0100 to 1200 UTC 30 Jul 2003. TheDA experiments differ in the variables assimilated into the model.Expt1\u201312-h forecastq (1.0 H11003 10H110024kg kgH110021)(N H11005 169)Bias MAE RMSE NRMSECTRL 18.883 24.465 29.723 1.0AT 17.605 24.051 28.801 0.969AQ 14.082 19.810 24.249 0.816ATQ 18.410 23.517 28.369 0.954TABLE 7. The summed total precipitation amounts (mm) overthe entire domain from experiments CTRL, TH1, TH2, TH3, andTH6.Expt 0\u20133h 3\u20136h 6\u20139h 9\u201312 hCTRL 96.9 0.0 0.9 7.6TH1 1095.6 0.0 0.5 5.0TH2 651.4 0.0 0.5 4.3TH3 653.4 0.0 0.5 4.3TH6 592.3 0.0 0.4 4.1MARCH 2007 DENG AND STULL 10476. Daily real-time runsThis section further tests MD_MT and the assimila-tion methodology for robustness over a longer periodby applying them in daily real-time forecast runs for 2months.a. Experimental settingThe real-time MC2 model at UBC has five one-wayself-nested grids with horizontal grid spacings of 108,36, 12, 4, and 2 km. The DA run is performed daily foronly the 2-km grid. The model time lines are shown inFig. 14. In a real-time setting, the grids at 108, 36, 12,and 4 km are actually run into the future at 1200 UTCon day 3. Limited by computational resources, the pri-mary real-time 2-km grid (hereafter the CTRL) drivenFIG. 12. Time series of observed surface temperature (dottedline with closed triangles) and the forecast temperature at thelowest model level from the control run (CTRL; dashed line withopen circles) and from experiment ATQ (solid line with closedcircles) for stations near McLure: (top) Kamloops ARPT,(middle) Sparks Lake, and (bottom) East Barriere.TABLE 9. The actual station elevations and the modeled stationelevations for the surface stations shown in Fig. 11.StationActualelev (m)Modeledelev (m)Elev error ofmodel (m)Kamloops ARPT 346.0 541.0 195.0Sparks Lake 1036.0 1030.0 H110026.0East Barriere 671.0 912.0 241.0Vernon 556.0 487.0 H1100269.0Kelowna College 300.0 399.0 99.0Penticton ARPT 344.0 569.0 225.0FIG. 13. Same as in Fig. 12, but for stations in Okanagan Valley:(top) Vernon, (middle) Kelowna College, and (bottom) PentictonARPT.1048 MONTHLY WEATHER REVIEW VOLUME 135by the 4-km output is run for only 27 h. This CTRL runis started at model-time 1200 UTC on day 1 (Fig. 14),and provides its 12-h output as a first guess for analysesat 0000 UTC on day 2 for the parallel DA run.The surface observations of potential temperatureand specific humidity at 0000 UTC on day 2 are ana-lyzed using the ADAS Bratseth (1986) scheme (Brew-ster 1996) modified to include the MD approach aftermountain-top refinement (MD_MT). The final analy-ses are obtained by combining the first guess, surface,and pseudo-upper-air data by using the PROF scheme.Out of the total 305 surface stations in the DA domain,22 stations were excluded from analysis and verificationbecause their actual and model elevation differed bymore than 500 m. Stencils of sharing factors at themodel grid points are generated for each of the remain-ing 283 surface stations using MD_MT, so that there isno need to recalculate the sharing factors duringeach day\u2019s analysis. These analyses are verified in sec-tion 6b.The 2-km DA parallel run is started at model time0000 UTC on day 2 and is integrated forward 15 h.Observation information are gradually incorporatedinto MC2 every 1200 s over a 1-h window using theIAU technique. Boundary conditions for both CTRLand DA runs were from MC2 4-km output. Verifica-tions of subsequent forecasts at the lowest model levelagainst surface observations are performed during the15-h forecast period for the DA runs, and are comparedwith verifications of the CTRL runs during the samevalid time. As discussed earlier in section 5, the CTRLruns are not strictly optimal. The improvement of theDA runs over the CTRL runs might be slightly over-estimated.Each day, the analysis obtains its first guess from the12-h output of the 2-km CTRL run, not from the 2-kmDA run. This implies that the local surface observationsassimilated are not incorporated into the forecast cycle;thus, each analysis and DA forecast are independent ofpast local surface observations. This daily cold start isgood for independent tests of the robustness of the DAmodule. Verification results for November and Decem-ber 2004 forecasts are given in section 6c. Due to ahardware failure of the supercomputer in December2004, the sample size for December is smaller than forNovember.b. Analysis resultsFigure 15 shows all 283 surface stations within theCanadian portion of the analysis and DA domain. Thenumber of observations varied from day to day depend-ing on the number of stations actually reporting. Be-cause verification was performed 1 day later than analy-sis (and some stations report late), the number of avail-able reports for verification in November 2004 waslarger than that for analysis (Fig. 16). December wassimilar (not shown).Five analysis methods are compared: GAUSS,TERR_DIFF, MD, MD_LSMG, and MD_MT. Formethod MD_LSMG, which included shoreline refine-ment (see DS05), 18 out of 283 surface weather stationsare identified to be stations over water, based on theland\u2013sea mask input data for the MC2 model at 2-kmgrid spacing.FIG. 14. Schematic diagram illustrating model time lines of thereal-time MC2 self-nested grids including the 2-km DA run. The2-km CTRL run is the regular real-time 2-km run. MC2 4-kmprovides boundary conditions for the 2-km CTRL run and theparallel 2-km assimilation run. SFC indicates surface data.FIG. 15. Surface weather stations (triangles) superposed on the2-km MC2 model topography (m). Darker shading indicateshigher elevations. Open triangles represent station locations withmissing temperature observations at this one sample analysis time0000 UTC 2 Nov 2004, while closed triangles indicate station lo-cations with available temperature observations at that analysistime.MARCH 2007 DENG AND STULL 10491) POTENTIAL TEMPERATURE ANALYSESTable 10 summarizes the verification statistics forNovember. The analyses from all methods have smallerpositive biases than the first guess (FSTG). The MAEand NRMSE of each method are largely decreasedfrom those of the first guess. Even though TERR_DIFF produces smaller bias than all other methods ex-cept MD_LSMG, it has the largest MAE and NRMSE.The MD approaches (MD, MD_MT, and MD_LSMG)have smaller MAE and NRMSE than GAUSS andTERR_DIFF, similar to the July 2003 case. TheMD_MT approach shows a minor increase of MAE andNRMSE. This is different from the July 2003 case.Comparatively, MD_LSMG has the smallest bias,MAE, and NRMSE, which is similar to what was foundfor the February 2003 case in DS05. But the improve-ment of MD_LSMG over MD is smaller, compared tothe February 2003 case. The verification statistics forDecember (not shown) show similar results to Novem-ber, but with the first guess having small negative biasand all analysis methods having smaller positive bias.2) SPECIFIC HUMIDITY ANALYSESThe bias, MAE, and NRMSE of the specific humidityanalyses and first guess for November are presented inTable 11. Each method gives largely reduced magni-tude of bias, MAE, and NRMSE compared to the firstguess. The MD approaches (MD, MD_MT, andMD_LSMG) produce better analyses with the smallernegative bias magnitude, MAE, and NRMSE thanGAUSS and TERR_DIFF, similar to the July 2003 case(Table 4). As there are no specific humidity observa-tions from the few mountain-top stations, MD_MT isequivalent to MD. MD_LSMG gives the smallest MAEand NRMSE of the specific humidity analyses. Decem-ber results are similar, but the magnitudes of bias,MAE, and NRMSE of FSTG are smaller than in No-vember. MD_LSMG slightly underperforms MD inDecember.c. Forecast results for MD_MT1) VERIFICATION OVER THE FULL FORECASTLENGTHFirst, look at the verification summary for differentnear-surface parameters over the full forecast lengthfrom the initial 0\u201315-h forecast period, over all surfaceTABLE 11. Same as in Table 10, but for verification of analyzedspecific humidity in Nov 2004.Method NBias (1.0 H1100310H110024kg kgH110021)MAE (1.0 H1100210H110024kg kgH110021) NRMSEFSTG 1125 H110025.987 7.800 1.0GAUSS 1125 H110020.451 1.892 0.292TERR_DIFF 1125 H110020.418 2.443 0.356MD 1125 H110020.169 1.687 0.270MD_MT 1125 H110020.169 1.687 0.270MD_LSMG 1125 H110020.464 1.553 0.247FIG. 16. Number of stations that had reports of (a) potentialtemperature and (b) specific humidity for each day of November2004. Each day in the x axis corresponds to \u201Cday 1\u201D in Fig. 14.Hence, the number of available reports are actually for a date thatis 1 day later than the date shown in x axis. Those indicated bytriangles are the numbers of reports used in analyses; whereasthose marked by diamonds are the numbers of reports used inverification.TABLE 10. Verification of analyzed potential temperatures inNovember 2004. For all these statistics, smaller magnitude isbetter.Method N Bias (K) MAE (K) NRMSEFSTG 4537 0.293 2.238 1.0GAUSS 4537 0.286 1.158 0.611TERR_DIFF 4537 0.261 1.278 0.644MD 4537 0.284 1.024 0.580MD_MT 4537 0.284 1.031 0.582MD_LSMG 4537 0.219 0.968 0.5631050 MONTHLY WEATHER REVIEW VOLUME 135stations that reported observations, and over the studyperiod of November and December 2004, respectively.Table 12 summarizes the verification statistics forNovember. For potential temperature forecasts, theDA run presents slightly larger negative bias than theCTRL run. The magnitude of the biases for both runs issmall. The MAE and NRMSE of the DA run are cor-respondingly smaller than those of the CTRL run. Themonthly averaged improvement in potential tempera-ture forecasts, as indicated by the value of NRMSE, issmaller compared to the July 2003 case study (see ex-periment ATQ in Table 6). A reduction in the improve-ment of near-surface potential temperature forecastsdue to assimilation of surface observations implies thatsurface thermal forcing in late fall may have beensmaller than in the summer.Humidity forecasts from the DA run have improvedbias, MAE, and NRMSE compared to the CTRL run.The improvement in humidity forecasts in November islarger when compared to experiment ATQ for the July2003 case study. The relatively dry July and wet No-vember help explain the difference in humidity forecastperformance.In contrast to the July 2003 case (see experimentATQ in Table 6), assimilating temperature and specifichumidity does not significantly improve the forecastquality of mean SLP for November.Assimilating temperature and specific humidity doesnot disturb the wind fields very much. This finding isdifferent from the July 2003 case study, perhaps be-cause performance is related to surface potential tem-perature improvement. When the correction to surfacepotential temperature is large (i.e., July 2003), the ini-tial imbalances between mass and wind fields due to asudden change in temperature fields are also large andcannot be removed by the IAU technique. However,when the correction to surface potential temperature issmall, the initial imbalances are also small and can besuppressed when the AIs are gradually introduced intothe model by the IAU technique.Verification statistics for December (Table 13) showthat potential temperature improvement is slightlylarger than for November, but is smaller compared toJuly 2003. The relative humidity forecasts from the DArun are less improved in December than in November.SLP and vector wind forecasts are neither improvednor degraded much in the DA run compared to theCTRL run.TABLE 13. Same as in Table 12, but for the month of Dec.Variable NBias MAE NRMSECTRL DA CTRL DA CTRL DAH9258 (K) 37 532 H110021.230 H110021.227 2.709 2.496 1.000 0.935RH(%) 21 198 2.685 3.127 10.612 10.272 1.000 0.983q (1.0 H11003 10H110024kg kgH110021) 10 056 H110024.053 H110023.581 7.436 6.413 1.000 0.888SLP (kPa) 10 878 H110020.055 H110020.055 0.336 0.334 1.000 1.002RMSVE NRMSVEVariable N CTRL DA CTRL DAV (m sH110021) 29 961 2.935 2.956 1.000 1.007TABLE 12. Verification statistics of near-surface parameters produced by the MC2 control (CTRL) and MD_MT DA runs. Thestatistics are calculated for all the observation\u2013forecast pairs over all stations that reported, and over the 15-h forecast period for themonth of Nov. The units for each variable are included for bias and MAE (or RMSVE). The NRMSE or NRMSVE is unitless.Variable NBias MAE NRMSECTRL DA CTRL DA CTRL DAH9258 (K) 55 230 H110020.426 H110020.524 2.286 2.160 1.000 0.959RH(%) 31 954 H110023.269 H110020.718 10.827 9.645 1.000 0.912q (1.0 H11003 10H110024kg kgH110021) 13 290 H110025.148 H110023.475 7.460 6.390 1.000 0.886SLP (kPa) 14 172 0.124 0.132 0.224 0.226 1.000 1.004RMSVE NRMSVEVariable N CTRL DA CTRL DAV (m sH110021) 43 823 2.696 2.693 1.000 0.999MARCH 2007 DENG AND STULL 10512) VERIFICATION BY FORECAST HOURTime series of monthly averaged NRMSE of surfacepotential temperature for November are shown in Fig.17a. An immediate improvement in NRMSE of the DArun is apparent at 1-h forecast, as new observation in-formation is introduced into the model within the 1-hwindow. Then, the improvement diminishes with fore-cast hour, but is still visible at 15 h. December (notshown) shows the same pattern as for November, butwith smaller NRMSE. The gradual decrease in the im-provement of the DA forecasts can be partly attributedto the effects of lateral boundary conditions advectinginto the domain, and partly to model errors that returnafter assimilating the observations.The DA improvement of specific humidity in bothNovember (see Fig. 17b) and December (not shown)behaves very similar to that of potential temperature.However, by assimilating temperature and specific hu-midity only, the model gives slightly degraded SLP andwind forecasts during the first several forecast hours forboth November (see Figs. 17c,d) and December (notshown). The deviation of the DA wind forecasts fromCTRL for November and December 2004 is muchsmaller than for July 2003. Near the middle and end ofthe forecast period, SLP and wind forecasts from theDA run are very similar to those from CTRL, againpossibly due to the assimilated information propagatingout of the domain, and due to mass and wind fieldsadjusted toward balance.7. SummaryThe MD method for analyzing complex terrain sur-face weather observations into a high-resolution NWPmodel (MC2) is tested with case studies and with 2months of daily real-time runs. The MC2 first guess,MD surface analysis, and pseudo-upper-air data (inter-polated from coarser-resolution analyses) are com-bined into a final 3D high-resolution analysis using aprofile method, which vertically spreads the surfacedata throughout the BL. The final analysis incrementsare incorporated into MC2 every 1200 s over a 1-h win-dow. This incremental analysis updating (Bloom et al.1996) exhibits better skill than instantaneous analysisincrement insertion.Also, a new mountain-top refinement (MD_MT)is proposed to the basic MD approach, and is testedalong with a previously suggested shoreline refine-ment (MD_LSMG). These approaches are also com-pared with Gaussian (GAUSS) and terrain differenceFIG. 17. Time series of monthly averaged NRMSE of surface parameter output from the CTRL and DA runsin November 2004: (a) potential temperature, (b) specific humidity, (c) mean SLP, and (d) vector wind.1052 MONTHLY WEATHER REVIEW VOLUME 135(TERR_DIFF) methods for surface data analysis in com-plex terrain. Verification against surface observations isdone both for the final analyses, and for the high-resolution forecasts started from MD_MT analyses.It is found that the MD approaches outperformGAUSS and TERR_DIFF when tested in steep moun-tainous regions of southern British Columbia, Canada.The mountain-top refinement adds increasing valueover MD as more mountain-top stations are included.The MD_LSMG shoreline refinement excels in coastalmountainous regions.Surface information assimilated at only the lowestmodel level is soon lost at the beginning of the forecastperiod. Larger improvement over the control (CTRL)run is achieved when surface information is spread up-ward throughout the BL. Combining the surface andpseudo-upper-air data gives slightly larger improve-ment over the CTRL forecast than assimilating onlysurface data throughout the whole BL. This impliesthat assimilation of surface data plays an important rolein reducing the model errors of near-surface weatherparameters, supporting the findings of Hacker and Sny-der (2005) who used an ensemble Kalman filtermethod.By assimilating surface temperature and specific hu-midity, the forecast quality of those parameters are bet-ter than CTRL. The DA improvement is the largest at1 h, and gradually decreases out to 15 h. The CTRL runin this paper is initialized by an analysis including olderobservations than the DA runs. This implies that theDA improvement over CTRL might be slightly over-estimated.The DA run gives poorer forecasts of near-surfacewinds and precipitation, which are not assimilated intothe model. This research addresses just a small part ofa much larger data assimilation need. Namely, it at-tempts to improve only the potential temperature andhumidity analyses and forecasts by assimilating thesedata. No attempt was made to improve precipitation orwind forecasts in the mountains. While there is indeeda critical need to improve forecasts of precipitation andwinds, those weather elements are much more difficultto assimilate successfully at high resolution in themountains, and no attempt is made to do it here.For these reasons, minimal discussion is given to pre-cipitation. The only reason precipitation is mentionedat all is so that the reader is aware what fields are andare not improved by the assimilation of only potentialtemperature and humidity.Future work could be done to assimilate sequences ofhourly surface observations in order to achieve betterquality for a longer forecast. Future investigationscould include the analysis of surface winds (probablythrough streamfunction and velocity potential) in com-plex terrain and the incorporation of winds, togetherwith temperature and specific humidity, into an NWPmodel. To get better precipitation forecasts while as-similating surface data, the assimilation method couldbe incorporated into a mesoscale three-dimensionalvariational data assimilation (3DVAR) system, whichcould also ease using an optimal CTRL.Acknowledgments. We thank J. Hacker from NCARfor his helpful recommendations and advice about thiswork, and S. Chamberland from RPN for MC2 usersupport. We also thank H. Modzelewski, G. Hicks II,and T. Cannon at UBC for their adept computer ad-ministration. Thanks also to two anonymous reviewers.Grant support came from the Canadian Natural Sci-ence and Engineering Research Council (NSERC), theBC Forest Investment Account, Environment Canada(EC &amp; LUTE), and the Canadian Foundation for Cli-mate and Atmospheric Science (CFCAS). GeophysicalDisaster CFD Center computers were funded by theCanadian Foundation for Innovation, the BC Knowl-edge Development Fund, and UBC.REFERENCESAlapaty, K., N. L. Seaman, D. S. Niyogi, and A. F. Hanna, 2001:Assimilating surface data to improve the accuracy of atmo-spheric boundary layer simulations. J. Appl. Meteor., 40,2068\u20132082.Barnes, S. L., 1964: A technique for maximizing details in numeri-cal weather map analysis. J. Appl. Meteor., 3, 396\u2013409.Benoit, R., M. Desgagne, P. Pellerin, S. Pellerin, Y. Chartier, andS. Desjardins, 1997: The Canadian MC2: A semi-Lagrangian,semi-implicit wideband atmospheric model suited for fine-scale process studies and simulation. Mon. Wea. Rev., 125,2382\u20132415.Bergeron, G., R. Laprise, D. Caya, A. Robert, M. Giguere, R.Benoit, and Y. Chartier, 1994: Formulation of MesoscaleCompressible Community (MC2) model. Internal Rep., Co-operative Centre for Research in Mesometeorology, 165 pp.Bloom, S. C., L. L. Takacs, A. M. da Silva, and D. Ledvina, 1996:Data assimilation using incremental analysis updates. Mon.Wea. Rev., 124, 1256\u20131271.Bratseth, A. M., 1986: Statistical interpolation by means of suc-cessive corrections. Tellus, 38A, 439\u2013447.Brewster, K. A., 1996: Application of a Bratseth analysis schemeincluding Doppler radar data. Preprints, 15th Conf. onWeather Forecasting and Analysis, Norfolk, VA, Amer. Me-teor. Soc., 92\u201395.\u2014\u2014, 2003: Phase-correcting data assimilation and application tostorm-scale numerical weather prediction. Part II: Applica-tion to a severe storm outbreak. Mon. Wea. Rev., 131, 493\u2013507.Deng, X., and R. Stull, 2005: A mesoscale analysis method forsurface potential temperature in mountainous and coastalterrain. Mon. Wea. Rev., 133, 389\u2013408.Hacker, J., and C. Snyder, 2005: Ensemble Kalman filter assimi-MARCH 2007 DENG AND STULL 1053lation of fixed screen-height observations in a parameterizedPBL. Mon. Wea. Rev., 133, 3260\u20133275.Laprise, R., D. Caya, G. Bergeron, and M. Gigu\u00E8re, 1997: Theformulation of the Andr\u00E9 Robert MC2 (mesoscale compress-ible community) model. Numerical Methods in Atmosphericand Oceanic Modelling: The Andr\u00E9 J. Robert Memorial Vol-ume, C. A. Lin, R. Laprise, and H. Ritchie, Eds., Atmo-sphere\u2013Ocean Series, Canadian Meteorological and Oceano-graphic Society/NRC Research Press, 195\u2013220.Miller, P. A., and S. G. Benjamin, 1992: A system for the hourlyassimilation of surface observations in mountainous and flatterrain. Mon. Wea. Rev., 120, 2342\u20132359.Ruggiero, F. H., K. D. Sashegyi, R. V. Madala, and S. Raman,1996: The use of surface observations in four-dimensionaldata assimilation using a mesoscale model. Mon. Wea. Rev.,124, 1018\u20131033.\u2014\u2014, G. D. Modica, and A. E. Lipton, 2000: Assimilation of sat-ellite imager data and surface observations to improve analy-sis of circulations forced by cloud shading contrasts. Mon.Wea. Rev., 128, 434\u2013448.Stauffer, D. R., N. L. Seaman, and F. S. Binkowski, 1991: Use offour-dimensional data assimilation in a limited-area meso-scale model. Part II: Effects of data assimilation within theplanetary boundary layer. Mon. Wea. Rev., 119, 734\u2013754.Stull, R. B., 1988: An Introduction to Boundary Layer Meteorol-ogy. Kluwer Academic, 666 pp.\u2014\u2014, 2000: Meteorology for Scientists and Engineers. 2d ed.Brooks/Cole, 502 pp.Xue, M., K. K. Droegemeier, and V. Wong, 2000: The AdvancedRegional Prediction System (ARPS)\u2014A multi-scale nonhy-drostatic atmospheric simulation and prediction model. PartI: Model dynamics and verification. Meteor. Atmos. Phys., 75,161\u2013193.\u2014\u2014, K. Brewster, D. Weber, K. W. Thomas, F. Kong, and E.Kemp, 2002: Real-time storm-scale forecast support forIHOP 2002 at CAPS. Preprints, 15th Conf. on NumericalWeather Prediction and 19th Conf. on Weather Analysis andForecasting, San Antonio, TX, Amer. Meteor. Soc., CD-ROM, 4B.3.Yee, S. Y. K., and A. J. Jackson, 1988: Blending of surface andrawinsonde data in mesoscale objective analysis. AFGLTech. Rep. 88-0144, Air Force Geophysics Laboratory,Hanscom AFB, MA, 31 pp. [NTIS ADA203984.]1054 MONTHLY WEATHER REVIEW VOLUME 135"@en
<https://open.library.ubc.ca/cIRcle/collections/facultyresearchandpublications/52383/items/1.0041839>