Form Approved REPORT DOCUMENTATION PAGE OMB No. 0104 0188 The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions lor reducing the burden, to the Department of Defense, Executive Services and Communications Directorate (0704-0188). Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number. PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ORGANIZATION. 3. DATES COVERED (From - To) 1. REPORT DATE (DD-MM-YYYY) 2. REPORT TYPE 24-11-2008 Journal Article 5a. CONTRACT NUMBER 4. TITLE AND SUBTITLE Value of Bulk Heat Flux Parameterizations for Ocean SST Prediction 5b. GRANT NUMBER 5c. PROGRAM ELEMENT NUMBER 0601153N 5d. PROJECT NUMBER 6. AUTHOR(S) Alan J. Wallcraft, A. Birol Kara, Harley E. Hurlburt, Eric P. Chassignet, George H. Halliwell 5e. TASK NUMBER 5f. WORK UNIT NUMBER 73-5732-16-5 8. PERFORMING ORGANIZATION 7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) REPORT NUMBER Naval Research Laboratory NRL/JA/7320-05-6069 Oceanography Division Stennis Space Center, MS 39529-5004 10. SPONSOR/MONITOR'S ACRONYM(S) 9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) ONR Office of Naval Research 800 N. Quincy St. Arlington, VA 22217-5660 11. SPONSOR/MONITOR'S REPORT NUMBER(S) 12. DISTRIBUTION/AVAILABILITY STATEMENT Approved for public release, distribution is unlimited. 13. SUPPLEMENTARY NOTES 14. ABSTRACT Hulk heat flux parameterization is an increasingly popular technique for forcing non-coupled ocean models. If sea surface temperature (SST) from the model is colder (warmer) than observed, then the net heat flux will be higher (lower) than observed; thus, bulk parameterizations tend to keep model SST close to observational SS I on long time scales. However, bulk parameterizations imply neither strong damping of SST variability nor strong relaxation to near-surface (e.g.. at 10 m) air temperature (Tj. This is demonstrated using SST simulations from a 0.72° * 0.72° cos(lat) (longitude * latitude) resolution global IIYhrid Coordinate Ocean Model (HYCOM) that does not include assimilation of any SST data or explicit relaxation to any SST climatology, but does use bulk heat fluxes. Results are discussed when climatological wind and thermal atmospheric forcing for HYCOM are constructed from three different archived numerical weather prediction (NWP) products: (I) the European Centre for Medium-Range Weather Forecasts (ECMWF) 15-year Re-Analysis during 197° -1993 (ERA-15). (2) ECMWF 40-year Re-AnaKsis 15. SUBJECT TERMS Bulk heat fluxes, ocean model SST, Exchange coefficients, atmospheric forcing, climate 17. LIMITATION OF 18. NUMBER 19a. NAME OF RESPONSIBLE PERSON 16. SECURITY CLASSIFICATION OF: ABSTRACT OF Alan Wallcraft a. REPORT b. ABSTRACT c. THIS PAGE PAGES UL 19b. TELEPHONE NUMBER /Include area code) Unclassified Unclassified Unclassified 18 228-688-4813 Standard Form 298 (Rev 8/98) Prescribed by ANSI Std. 239.18 • Available online at www.sciencedirect.com JOURNAL OF ScienceDirect MARINE SYSTEMS ELSEVIER Journal of Marine Systems 74 (2008) 241 258 - www.elsevier.com/locate/jmarsys Value of bulk heat flux parameterizations for ocean SST prediction 3 a Alan J. Wallcraft"*, A. Birol Kara , Harley E. Hurlburt , b c Eric P. Chassignet , George H. Halliwell " Oceanography Division. Naval Research Laboratory, Stennis Space Center, MS, I SI ('enter for Ot can Atmospheric Prediction Studies, ami Department of Oceanography, Florida State University, Tallahassee. Florida. USA L Rosenstiel School of Marine and Atmospheric Science. University of Miami, FL. i'SA Received 30 March 201)7; received in revised form 4 January 200X: accepted 11 January 2008 Available online 10 March 2008 Abstract Bulk heat flux parameterization is an increasingly popular technique for forcing non-coupled ocean models. If sea surface temperature (SST) from the model is colder (wanner) than observed, then the net heal flux will be higher (lower) than observed; thus, bulk parameterizations tend to keep model SST close to observational SST on long time scales. However, bulk parameterizations imply neither strong damping of SST variability nor strong relaxation to near-surface (e.g., at 10 m) air temperature (T,). This is demonstrated using SST simulations from a 0.72° x 0.72° cos(lat) (longitude * latitude) resolution global HYbrid Coordinate Ocean Model (HYCOM) that does not include assimilation of any SST data or explicit relaxation to any SST climatology, but does use bulk heat Muxes. Results are discussed when climatological wind and thermal atmospheric forcing for HYCOM are constructed from three different archived numerical weather prediction (NWT) products: (1) the European Centre for Medium-Range Weather Forecasts (ECMWF) 15-year Re-Analysis during 1079-1993 (ERA-15), (2) ECMWF 40-year Re-Analysis (ERA-40) during 1979-2002, and (3) the Common Ocean Reference Experiment Corrected Normal Year forcing version 1.0 (CORE-CNY) based on the National (enters for Environmental Prediction (NCEP) re-analysis which spans 1948-2002. To investigate the implications of the bulk heat flux approach as relaxation to SST and 7" , HYCOM SST simulations are used to demonstrate that model SST errors with respect to the a National Oceanic Atmospheric Administration (NOAA) SST climatology do not look like T„ SST fields from NWP products Such a result is confirmed for all simulations forced with ERA-15, ERA ^40 and CORE-CNY, separately. Overall, global averages of mean HYCOM SST error are 0.2 °C (1.5 °C), 0.4 °C (1.7 °C) and 0.6 °C (2.3 °C) with respect to NOAA SST (NWP T ) climatology when a the model uses atmospheric forcing from ERA-15, ERA-40 and CORE-CNY, respectively. C 2008 Elsevier B.V. All rights reserved. Keywords: Hulk heal fluxes; Ocean model SST; Exchange coefficients; Atmospheric forcing: Climate I. Introduction and motivation * Corresponding author , fi -OTa;7a</</n?.w.v.alan.wallcraft^nrlssc.navy.mil(A.J.Wallcraft), Ocean general circulation model (OGCM) simula- bm.i.kara.,/nrissc.naNy.mil (A.B.Kara). g || performed using prescribed atmo- tjons are encra y harlcy.hurlburtff/.nrlssc.navy.mil (II.E. llurlburt), .._._.. . _ . ' , . .,.„,.. ,, spheric forcing fields, namelv, momentum tliix (c.o., r e J e cchassigiietW;eoaps. tsu.edu (F..P. C nassignet), v o» ghalhwdl«, rsmas miami.edu (G ll Halliwell). wind stress) and scalar forcing (e.g., net shortwave and URL: http:/'www732o.nrissc navy.mil (A J Wallcraft). longwave radiation at the sea surface, air temperature 0924-7963/$ - see front matter C 2008 Elsevier B.V All rights res' doi:10.IOI6/j.jmarsys.2008.0l.009 20090306220 242 A.I. WaUcrafi el al. Journal of Marine Systems ^4 (2008) 241 25S and air mixing ratio at 10 m above the sea surface). the exchange coefficients for latent and sensible heat These forcing fields are typically obtained from fluxes (e.g., Kara et al.. 2002). The most important criticisms of bulk parameteriza- archived NWP products. Examples of commonly-used NWP products include the European Centre for tions are that (1) an infinite heat capacity exists between Medium-Range Weather forecasts (ECMWF) 15-year ocean and atmosphere, i.e., it does not give a closed Re-Analysis (Gibson et al., 1999), ECMWF system, (2) since the bulk formula can be linearized as a TERA-15) 40-year Re-Analysis (ERA-40) (Kallberg et al.. 2004), function of a difference between equilibrium tempera- National Centers for Environmental Prediction (NCEP) ture and SST (Haney, 1971; Han. 1984: Paiva and Re-Analysis (Kanamitsu el al., 2002), and the Fleet Chassignet, 2001), it can smear the model SST when Numerical Meteorology and Oceanography Center the model resolution is finer than the forcing, and (FNMOC) Navy Operational Global Atmospheric (3) The difference, i.e., T SST. generally represent a Prediction System (NOGAPS) (Rosmond et al., 2002). the boundary layer physics of the NWP products, but All NWPs mentioned above use relatively sophisti- may be a poor approximation to observations. Note cated boundary layer sub-models to calculate sur- that the statement in (1) typically holds if the atmo- face fluxes, but they still depend on SST which is sphere temperature is kept constant over a longer period, usually an analyzed (non-prognostic) field. The SST for example if climatological values arc used. Otherwise used is now accurate enough that errors in surface fluxes the atmosphere temperature exhibits a time tendency too, which should reflect the effect of the heat flux. at the NWP grid-scale are dominated by errors in atmospheric fields, such as wind speed and cloudiness. In this paper, we examine the impact of bulk heat (lux These errors can be large, and even the best surface parameterization on the climatological SST produced by heat flux fields from all sources (NWP products or simulations from a particular global OGCM, with no observation-based climatologies) do not give a closed assimilation of, or explicit relaxation to. any observed SST climatology. In particular, reasons for preferring the global heat budget. In addition, using climatological estimates of total heat flux to force an ocean model bulk parameterization rather the net heat flux plus an usually results in unrealistic model SST, presenting an explicit relaxation to SST arc discussed. inconsistency with the imposed surface flux (Bamier A few salient reasons for using SST from an OGCM et al., 1995). to discuss the bulk heat flux versus the direct relaxa- Even perfect fluxes, with a closed global heat budget, tion approach are (I) SST plays an important role in cannot be used as the only forcing for a stand-alone air-sea interactions through the net heat flux at the ocean model because the model's "climatology" is not sea surface (e.g., Alexander and Scott, 1997); (2) it is perfect, so such fluxes will lead to SST drift (e.g., typically used in the bulk heat flux formulation to Hughes and Weaver, 1996). The conventional approach parameterize the stability (Kara et al., 2000); (3) it to correcting this drift is to add relaxation to observed plays a major role in controlling long term climate SST, either in place of, or in addition to, prescribed total variations, such as the North Atlantic Oscillation heat (luxes (e.g., Bamier et al., 1995; Maltmd et al., (Pacth et al., 2003); (4) it is one of the best observed 1998). As stated by Killworth et al. (2000), researchers variables of the upper ocean over the global ocean; therefore, SSTs simulated from an OGCM can be at this time did not possess consistent surface fields with which to force their models. However, the use of a bulk easily validated over the globe (Schopf and l.oughe. 1995). In addition, SST from an atmospherically- parameterization has become an alternative approach because it requires no explicit relaxation to observed forced OGCM (no assimilation of any ocean data, oceanic quantities (e.g., Wallcraft et al.. 2003; Kara including SST) is used because one of our goals is to el al., 2003). In addition, using bulk parameterizations demonstrate that, although the bulk heat flux formula- (e.g., Parkinson and Washington. 1979), regional tion forcing an OGCM surface temperature field coupled ice-ocean models, such as those in the Baltic includes the knowledge of air temperature near the Sea, have been integrated over decades (Lchmann and sea surface, such use docs not imply an improper form Hinrichscn. 2000). of SST restoring within the model. Bulk parameterizations use the difference 7 SST, The paper is organized as follows. Section 2 a where 7[, is the air temperature and SST the ocean model discusses the traditional bulk heat flux approach for sea surface temperature. The air temperature and other OGCMs. Section 3 gives details of the OGCM used in parameters (air density at the air-sea interface, mixing this study. Section 4 presents SST results from an ratio, and wind speed at 10 in above the sea surface) are OGCM in relation to the bulk heat flux parameteriza- typically from NWP products and are used to calculate tion, and Section 5 gives conclusions of this paper. ,/../. Walkrafi et al. I Journal ofMarine Systems 74(2008)241 25/1 !43 coefficients {C, and C ) used in HYCOM include the 2. Bulk heat flux approach s effects of boundary layer stability, and are expressed as The total heat flux (Q ) at the sea surface can be polynomial functions of T SST, v.„ and q q , the ncl a a s expressed as follows: latter through relative humidity (RH), where RH is at the air-sea interface (see http://www7320.nrlssc.navy.mil/ nasec/nasec.html). In the polynomial functions, the (I) (?nc. = C?sw " 01* + Ql + Qs, effects of water vapor flux in calculating the exchange where (? is the net shortwave radiation at the sea surface, coefficients are taken into account through RH effects svv that are especially important at low wind speeds (Kara Qi is the net longwave radiation at the sea surface, Q K t is the latent heat flux, and Q is the sensible heat flux. et al., 2005d). The total heat flux (Eq. (1)) is therefore s All these individual components, and their total, Q „ expressed as a polynomial function of the model SST nc arc typically available from the archived real-time and re- and can be linearized to first order as analysis data sets or from climatologies (ERA-15, ERA-40 Q, - ;.(T* SST). and CORE-CNY are used here). Shortwave, (9 and (4) Kt SWS longwave. Q^ radiation can be directly measured at the surface, but there are far too few such observations to in which the apparent equilibrium temperature T* and constrain even regional fluxes. Thus, in practice all the relaxation coefficient '/. are time and space components arc estimated. All components are affected dependent and can be computed from the atmospheric by surface type i.e., land versus sea versus sea-ice (Gamitt fields (Haney, 1971: Han. 1984; Paiva and Chassignet, ct al., 1998), but here wc will only consider the open ocean. 2001). We have included Q in the Haney parameter- sw Both shortwave and longwave radiation depend heavily on ization, as is usually done, even though it does not cloudiness, and all flux components except shortwave depend on SST. The total heat flux therefore implies a radiation depend on SST (e.g., (iill 1982: Gleckler and temperature relaxation, but not toward the atmospheric Weare, 1997). HYCOM uses a penetrating solar radiation temperature r , nor toward a "correct" SST (7^). The a scheme (Kara ct al., 2005a) that accounts for spatial and equilibrium temperature T* indeed must differ from T , c temporal water turbidity (Kara et al.. 2005b,c). as it would otherwise imply zero heat flux when the While SST drift in OGCMs may be significantly model SST is equal to T . c reduced by relaxation to observed SST (e.g., Barnier et al., 1905), such an approach has the important short- 2.1. Preference for a bulk parameterization coming that the correct result is prescribed in the model simulation. Rather than using a direct SST relaxation Instead of using a direct relaxation to any SST scheme, an alternative approach is to calculate fluxes climatology, a common alternative has been to use N WP from a bulk parameterization based on near-surface net flux plus an explicit relaxation to the correct SST, 7" . c atmospheric fields using only near-surface wind speed Since SST is one of the best observed geophysical fields, (v.,). or also including air mixing ratio (i/. ), near-surface why should a bulk parameterization be preferred over t air temperature ('/],) (all of which arc 10 m above the sea the NWP net flux plus an explicit relaxation to SST? surface), mixing ratio for sea water (</ ) and model SST, One answer is that the c-folding time implicit in accurate s as explained in Kara et al. (2005d). bulk parameterizations is representative of the air-sea Bulk paramctcrizations arc based on statistical fits to interface. It is possible to construct explicit relaxation observations. They arc usually inexpensive to calculate, terms that are patterned on bulk paramctcrizations. or and the best paramctcrizations arc generally very alternatively to use measured climatological e-folding accurate. Examples of typical bulk paramctcrizations times from observations or NWP fields (or the results of for latent and sensible heat fluxes can be found in a bulk parameterization applied to NWP fields). These DcCosmo ct al. (1996). Kara ct al. (2000, 2002). and approaches are all of the Haney type (Haney. 1971; I Ian, 1984; da Silva et al.. 1994; Paiva and Chassignet. 2001) Fairall ct al. (2003). The formulation used in HYCOM is briefly described as follows: as follows C?s = c\r > v (7; -SST) (2) SST). £(SST) = g(T )+^| .(r p/ a il 1 s n c (3) Ql - C'//./' V ,((/a </J M ; where the effective SST relaxation e-folding time scale implied by 'ff \ depends on many factors. However, it The air density at the air-sea interface (p ) in kg m r a is certainly positive and typically large enough that is determined using the ideal gas law. The exchange 244 I.I. Wallcrafl el al. IJournal ofMarine Systems 74(2008)241 258 model SST docs not exhibit long term drift. In practice the SST, therefore introducing oceanic information into relaxation terms are often simpler than this, i.e., a the heat flux formulation, the long term (climatological) constant c-folding time and some approximation to T . mean of the difference between near-surface 7" and SST c a from an observation-based data set and the SST used by Eddy-resolving, atmospherically-forced ocean models without ocean data assimilation should not relax strongly NWP products arc investigated over the global ocean to synoptic observed SST because the observed SST will (Fig. I). Specifically, the observation-based data set is include SST anomalies associated with eddies and the Comprehensive Ocean Atmosphere Data Set western boundary currents at different locations than in Climatological mean bias: T<, 1\ the model simulation. This is because of flow instabilities in the model which are not a deterministic response to the forcing. The bulk parameterization does not use the NWP SST (7' ) explicitly, and the NWP fields arc often on a s much coarser grid than the ocean model one that barely resolves such details as eddy and ocean current locations. However, this raises the question of to what degree the implied relaxation in the heat flux can smear the model SST when the model resolution is finer than the forcing. It also raises the question as to whether the NWP atmo- spheric temperature /,, is strongly correlated to the NWP SST 7" , therefore introducing oceanic information in the s heat flux formulation. In contrast, eddy-resolving stand- alone ocean models that assimilate sea surface height (from satellite altimeters) have ocean fronts and eddies in the observed location (Smedstad ct al., 2003; Chassignet et al., 2006; Shriver et al., 2007) and therefore can relax directly to high resolution observed SST, although they may instead (or in addition) directly assimilate SST. From the preceding discussion, simple relaxation is never appropriate as the only heat flux term, so it must be augmented either with the NWP net heat flux or a bulk parameterization heat flux. Regions where the ocean model's mean SSTs arc less accurate can often be traced to deficiencies in the atmospheric fields. In particular, systematic biases in wind speed will always lead to poor fluxes from the bulk parameterization and biases in cloudiness will have a large effect on shortwave radiation (partially offset by changes to downward longwave radiation). Biased winds are unlikely to provide better surface heat fluxes from the sophisticated boundary layer sub-model in the NWP than from the bulk parameterization. Assimilation of wind speeds from satellites should improve the accuracy of NWP surface winds, but there are still large biases in some coastal regions, and there may be smaller systematic biases on larger scales. Cloud cover is difficult to obtain accurately, and known to be deficient regionally in both archived real-time and re-analysis products. Fig. I. Climatological mean of the difference between near-surface air 2.2. Relationship between air temperature and SST temperature (7" ) and sea surface temperature (i.e., T,- I'J over the a global ocean Both T, and 7; fields from COADS, FRA-15. ERA-40 In order to answer the question as to whether the and CORE-CNY are interpolated to the 0.72" IIYCOM grid, and the atmospheric temperature, ?",, is strongly correlated to difference between the two is formed at each grid point A..I. Walkraft et al. /Journal of Marine Systems 74 (200X1 241 258 245 (C'OADS) climatology based largely on ship observa- corrected for excessively cold values in the Antarctic. Ignoring high latitudes where sea-ice forms, the fields tions and buoy measurements during 1945-1989 (da Silva et al.. 1994) and the NWP data products are the re- are broadly similar to each other with climatological T s wanner than 7" nearly everywhere but usually by analysis products, ERA-15 (1979-1993), ERA-40 a (1979-2002) and CORE-CNY (1948-2002). COADS < 1 °C. The exception is in areas where strong currents is the only observation-based climatology among these transport water which is significantly warmer than the products, and is intended to complement the compar- air in the mean, e.g.. Gulf Stream, Kuroshio and Agulhas. Further details can be found in Kara et al. isons. This is the new 1/2° x 1/2° climatology based on (2007). the atlas of surface marine data. Supplement B (http:// The same long term mean fields shown in Fig. 1 arc- www.nodc.noaa.gov/OC5 'asmdnew.html). Note that hereinafter the COADS SST, as well as the SSTs used plotted as scatter diagrams of 7/, versus T over the s in the NWP products will be denoted as 7*. global ocean (Fig. 2). The regions where ice is located Although the time periods over which each climatol- are excluded in the evaluation procedure. There is a strong linear relationship between T., and SST with R ogy (COADS, ERA-15. ERA-40 and NCEP) is constructed are different, they arc all > 15 years which values >0.99 for all data sets. This is clearly evident is adequate to represent a long term mean over the from the slope values being close to I in all cases (1.008, global ocean. Near-surface T„ and /' data for each 0.991, 1.012 and 1.033 for COADS, ERA-15, ERA-40 s product are available from the National Center Atmo- and CORE-CNY, respectively). Intercept values arc also spheric Research (NCAR) web site (available online at similar with values of 0.468, 0.911, 0.711 and 0.405. http://dss.ucar.edu,catalogs). COADS provides The difference between 7/, and T (i.e., 7"., T ) ranges s s from - 0.64 °C for COADS to - 1.09 °C for CORE-CNY monthly mean 7" - 7" fields directly, while for the u s (Table I). three NWP products we form climatological means using T and the T fields archived at sub-daily intervals. In contrast, near-surface T and 7" 7" arc only a s a a s Near-surface atmospheric temperature 7~, from ERA-15, weakly correlated (Fig. 3) with linear correlation coef- ERA-40 and NCEP is at 2 m, while that from COADS is ficients of 0.12, 0.14. 0.17 and 0.25 for COADS, at 10 m above the sea surface. The adjustment from 2 to ERA-15, ERA-40 and CORE-CNY. None of the cor- 10 m T is very small as indicated by the comparison of relation values arc statistically significant in comparison A 7], values in fable I and hence was ignored. to a correlation value of 0 at a 95% confidence interval When sea-ice is present, surface T is over sea-ice in based on the student's Mcst. Thus, / itself does not a ;1 the re-analysis products, but not in COADS. The 7j, field have strong influence on T 7" . a s from NCEP/NCAR is actually obtained from the CORE-CNY (Corrected Normal Year) data set which 3. HYCOM description The HYbrid Coordinate Ocean Model (HYCOM) is rable 1 based on a primitive-equation formulation discussed by Global averages tor near-surface air temperature ( /',). T, and difference between the two from the observational-based C'OADS data set and Bleck (2002), in detail. There have been several three NWP-based products HYCOM applications investigating a variety of pro- cesses in different ocean basins and enclosed seas. T, 7", Data 7, 7; Examples of such studies include the North Atlantic 22.0 21.4 COADS 0.6 (Chassignel el al.. 2003; Halliwell 2004; Thackcr et al., 21.4 22 1 ERA-15 -0.7 F.RA-40 21.2 22.1 -0 9 2004), the Indian Ocean (Han et al.. 2004). the tropical COARE-CNY 20.9 22.0 1.1 Pacific (Shaji et al.. 2005), and the Black Sea (Kara Further information about COADS, ERA-15, ERA-40 and CORE- et al.. 2005a,b,c). HYCOM has also been used as the CNY can be found online at several web sites. For example, as of this ocean component of a coupled atmosphere ocean writing, the web addresses are http: iridlldeo.columbia.edu/ model in global climate studies (e.g.. Sun and llansen SOURCES/.DAS1LVA SMD94 for COADS. http://badc.nere.ac.uk,' 2003: Bleck and Sun 2004), and for cddy-frcsolving data'ccmwf-era'ERA.htitil for HRA-15. http://badc.nerc.ac.uk/data ocean prediction (Chassignet et al.. 2006). ecmw'f-e40c40.. background html for HRA-40, and http: datal.gfdl. noaa.gov/nomads/forms.mom4 (OKI. CNYF IpO.html for CORE- For the present study, we will introduce FIYCOM CNY. The web addresses may be subject to changes by the originators. configured for the global ocean, including the latest All the NWP products have different boundary layer paramctcriza- advances in the model development (sec Appendix for lions. physics, data assimilation methods and different satellite data details). The model presented here is a stand-alone used m the assimilations. Therefore, differences in their output ocean model with no assimilation of any ocean data, variables arc expected. A.J. Wallcraft el al. Journal of Marine Systems 74 (2008) 241 258 246 including SST, and no relaxation to any other data except sea surface salinity (SSS) to keep the evaporation minus precipitation balance on track. General features of the global atmospherically-forced HYCOM are given in 4 ' 1 ' 1 1 \s 30 _ COADS * — 20 — 10 \i2^r 0 , -h 7T i 1 i I 10 20 30 4 ' 1 ' 1 \i . L/ 30 ^^^F' - ERA-15 — 20 G ' - — ^' 10 , t- , 1 , 1 0 TK 10 20 30 0 4 1 ' 1 30 - ERA r _ 2n c — — - - '•—• — — ^ 10 1 , 0 f- 7r 1 4 10 20 30 2 jj 4 ' 1 ' 1 l [£ 30 CORE-CNY 0 1 h, " — 20 j/KF G 3 -2 s- - '— — ^' 10 "W 0 10 20 30 if, 1 , 1 , +" 0 7-,»C, 20 30 0 T ( C) a fig. 3. Scatter plots of near-surface air temperature (7;,) versus f, - 7", based on annual mean values (see fig I), each at I ° bins over the Fig. 2 Scatter plots of T., versus r„ based on annual mean values, each global ocean. at 1° bins over the global ocean. Both 7], and T, fields are from COADS, KRA-15. ERA-40 and CORK-CNY, respectively. A..I. Wallcraft el al. I Journal ofMarine Systems 74 (2008) 241 25S 247 Section 3.1 and the atmospheric forcing used in the density with depth between isopycnal coordinate model simulations is explained in Section 3.2. surfaces are based on the 1/4° Generalized Digital Environmental Model (GDEM) climatology (NAVO- 3.1. General features of global HYCOM CEANO, 2003). The density difference values were chosen so that the layers tend to become thicker with IIYCOM contains five prognostic equations: two for increasing depth, with the lowest abyssal layer being the the horizontal velocity components, a mass continuity or thickest. The near-surface z-levcl regime is a natural layer thickness tendency equation and two conservation consequence of HYCOM's minimum layer thickness. In equations for a pair of thermodynamic variables, such as this case, the minimum thickness of layer 1 is 3 in, and salt and potential temperature or salt and potential this minimum increases 1.125* per layer up to a density (Bleck, 2002). The model behaves like a maximum at 12 m, and target densities arc chosen, conventional o (terrain-following) model in very such that at least the top four layers are in r-level shallow oceanic regions, like a z-level (fixed-depth) coordinates. coordinate model in the mixed layer or other unstratified The simulations use realistic bottom topography regions, and like an isopyenic-coordinate model in constructed from the NRL 2 min resolution bathymetry. stratified regions. However, the model is not limited to Extensive quality control of bottom topography was these coordinate types (Chassignet et al., 2003). performed in straits and near coastlines. Monthly mean Typically, HYCOM has isopycnal coordinates in the temperature and salinity from the GDEM climatology in stratified ocean but uses the layered continuity equation August are used to initialize the model. There is a to make a dynamically smooth transition to r-lcvcls in relaxation to monthly mean SSS from the Polar science the unstratified surface mixed layer and rr-lcvels in center Hydrographic Climatology (PHC). The PHC shallow water. The optimal coordinate is chosen every climatology is chosen for its accuracy in the Arctic time step using a hybrid coordinate generator. Thus, the region (Stcclc et al., 2001). The reference mixed layer model automatically generates the lighter isopycnal thickness for the SSS relaxation is 30 m (30 days in 30 in layers needed for the pycnocline during summer, but c-folding time). The actual e-folding time depends on during winter the same layers define fixed-depth z-level the mixed layer depth (MLD) and is 30 * 30/MLD days, coordinates. i.e. it is more rapid when the MLD is shallow and less HYCOM domain used in this paper spans the global so when it is deep. Such a relaxation is necessary to ocean from 78°S to 90°N. It has a 0.72° equatorial prevent SSS drift, and is in addition to the evaporation Mercator grid between 78°S-47°N, with an Arctic bi- precipitation budget. polar patch above 47°N. Average zonal (longitudinal) HYCOM has several mixed layer/vertical mixing grid resolution for the 0.72° global model varies from options (see Halliwell (2004) for a discussion and w 80 km at the equator to « 60 km at mid-latitudes (e.g., evaluation). In this paper, the non-slab K-Profile at 40° N). The meridional (latitudinal) grid resolution is Parameterization (KPP) of Large et al. (1997) is used. doubled near the equator to better resolve the equatorial wave-quide and halved in the Antarctic for computa- 3.2. Atmospheric forcing tional efficiency. Hereinafter, the model resolution will be referred to as 0.72° for simplicity. Zonal and The model reads in the following time-varying meridional array lengths are 500 and 457. respectively. atmospheric forcing fields: wind forcing (zonal and At this resolution, coastal regions are not represented in meridional components of wind stress, wind speed at great detail (so sigma-levels are not used), and the model 10 m above the sea surface) and thermal forcing (air land sea boundary is at the 50 m isobath (with a closed temperature and air mixing ratio at 10 m above the sea Bering Strait). surface, precipitation, net shortwave radiation and net All global simulations use sigmaO (sigma-theta), i.e. longwave radiation at the sea surface). The wind/ potential density referenced to the surface (0 dbar) with thermal forcing was constructed from three NWP no thennobaric correction. There are 26 hybrid layers in products: (i) 1.125°* 1.125° ERA-15, (ii) 1.125°* the vertical. In general, the model needs fewer vertical 1.125° ERA-40 and (iii) 1.875° x 1.875° CORE-CNY. coordinate surfaces than, say, a conventional r-lcvel When using each product, the goal is to include high model, because isopycnals arc more efficient in temporal variability in the forcing, while maintaining a representing the stratified ocean, as discussed in Hurlburt climatology. et al. (1996) and Kara et al. (2005a). The target density The original atmospheric forcing data set from ERA-15 values for the isopycnals and the decreasing change in spans the period I January 1979-31 December 1993. A.J. Wallcraft el al. i Journal of Marine Systems 74 COOS) 241 258 248 The climatologic.il ERA-15 forcing used here consists of coefficient at 490 nm) data set from Sea-viewing Wide 15-year monthly averages over 1979-1993, interpolated in Field-of-view Sensor (SeaWiFS) during 1997-2001. time to 6 hourly and with 6 hourly sub-monthly anomalies Thus, using ocean color data, the effects of water turbidity from operational ECMWF in September 1994 to Septem- are included in the model simulations through the ber 1995 added to the winds only. Choosing another time attenuation depth (A> in m) for shortwave radiation. AR period for the 6 hourly wind anomalies (other than 1994 The rate of heating/cooling of model layers in the upper 95) did not make any significant impact on the model SST. ocean is obtained from the net heat flux absorbed from the HCMWF used a spectral model to generate the ERA-15 sea surface down to a depth, including water turbidity dataset. The ERA-15 re-analysis project incorporated a effects (e.g., Kara ct al., 2005a). Previously, it was shown number of in-situ and satellite-based data types over the 15- that water turbidity can be quite significant in SST year period. SST analyses for ERA-15 are provided by the simulations, especially in equatorial regions (e.g., Kara United Kingdom Meteorological Office (UKMO) for the ct al., 2004. 2005b). early period, and by NOAA from November 1981 The net surface heat flux that has been absorbed (or onwards. Sea-ice cover has been derived from satellite data. lost) by the upper ocean to depth is parameterized as the The ERA-40 project applies a modem variational sum of the downward surface solar radiance, upward data assimilation technique for the past conventional longwave radiation, and the downward latent and and satellite observations. The model physics and the sensible heat fluxes (see Section 2). Net solar radiation surface parameterization have been upgraded and (the sum of net shortwave and longwave radiation) at the improved since ERA-15. A significant difference sea surface is so dependent on cloudiness that it is taken between ERA-40 and ERA-15 is in the use of satellite directly from the given NWP product (ECMWF or data. ERA-40 uses the Advanced Tiros Operational CORE-CNY) for use in the HYCOM. The net longwave Vertical Sounde (ATOVS) radiances directly, while in flux is the sum of downward longwave (from the ERA-15 temperature and humidity retrievals were used. atmosphere) and upward black-body radiation. The In addition. Special Sensor Microwave Imager (SSM/1) NWP (input) biack-body radiation is corrected within and European Remote Sensing Satellite (ERS) data are HYCOM to allow for the difference between NWP SST, used in ERA-40. The SST/Ice data set produced by the 7" , and HYCOM SST (Kara ct al., 2005b). The s Hadley Centre and National Oceanic and Atmospheric downward longwave radiation is often not archived, Administration (NOAA)/National Environmental Satel- but if archived, input net longwave is calculated using lite. Data, and Information Service (NESD1S) has been an archived NWP SST made available to the ERA-40 project. For use in Latent and sensible heat fluxes at the air sea interface the analyses and HYCOM simulations performed here, are calculated using computationally efficient bulk 25-year monthly climatologies of atmospheric forcing formulae that include the effects of dynamic stability parameters were formed from the ERA-40 archives (Kara et al.. 2005d). The details of the parameterizations spanning 1979-2002 with the same 6 hourly wind are given in Section 2. HYCOM treats rivers as a "runoff" anomalies added as in the ERA-15 case. addition to the surface precipitation field. The flow is first The atmospheric forcing from CORE-CNY is for a applied to a single ocean grid point and smoothed over single year and combines re-analysis data from NCEP surrounding ocean grid points, yielding a contribution to with satellite measurements to reduce errors existing in precipitation in m s~ . the original NCEP fields, especially radiation fields (e.g., Lee ct al., 2005). Heat fluxes for the CORE-CNY 4. HYCOM SST simulations forcing are produced by applying the NCAR bulk parameterization to NOAA SSTs, however HYCOM In this section, we investigate whether or not monthly calculates its own latent and sensible heat fluxes based mean SSTs obtained from atmospherically-forced on model SST (see below). The forcing includes 6 HYCOM, which includes the effects of bulk parameter- hourly representative variability in all its fields. izations in its surface energy balance (see Section 2), is In addition to wind and thermal forcing, HYCOM strongly controlled by the near-surface T used in the a forcing incorporates monthly mean climatologies of river atmospheric forcing. The global HYCOM simulations discharge and satellite-based attenuation coefficients for presented in this paper were performed with no assimila- Photosynthetieally Active Radiation (A> in m '). The tion of any oceanic data except initialization from AR shortwave radiation at depth is calculated using a spatially climatology. There is only weak relaxation to sea surface and temporally varying monthly Ap climatology as salinity to keep the evaporation-precipitation budget on AK processed from the daily-averaged k (attenuation track in the model. Model simulations arc performed 490 A.l. Walkraft el at. /Journal of Marine Systems 74 (2008) 241 2SH 249 using atmospheric wind/thermal forcing from each NWP arc mainly designed for large-scale climate studies. product (ERA-15, ERA-40 and CORE-CNY), separately. They arc derived by a linear interpolation of the weekly We did not perform any simulations using the alternative optimal interpolation (01), and use in-situ and satellite approach of applying NWP net heat fluxes and relaxing to SSTs, making it a reliable candidate for model-data "observed" SS1, both because this explicitly includes the comparisons over the global ocean. The existence of the "answer", at least for SST, in the model run and because ice field in the NOAA data set is also an advantage for we would have to choose one particular relaxation scheme the model SST validation in the Arctic and Antarctic. which would still leave open the question of whether We will also use 7" climatologies from NWP products to a some other variant would be more appropriate. compare with HYCOM SSTs in Section 4.2. Near statistical equilibrium was reached for each Monthly mean HYCOM SST obtained from simulation after four model years. A linear regression HYCOM simulations using wind and thermal forcing analysis was performed for domain averaged quantities from ERA-15, ERA-40 and CORE-CNY is validated (layer temperature, salinity, potential and kinetic energy, using mean error (ME), root-mean-square (RMS) SST etc.) to investigate statistical equilibrium in each layer, difference, correlation coefficient (R) and non-dimen- and is expressed numerically as % change per decade. sional skill score (SS). Let A', (/^ 1, 2, , 12) be the set of The model is deemed to be in statistical equilibrium monthly mean NOAA reference (observed) SST values when the rate of potential energy change is acceptably from January to December, and let Y, (i- I, 2, , 12) be small (e.g., < 1% in 5 years) in all layers. The statistical the set of corresponding HYCOM estimates at a model equilibrium was accomplished using climatological grid point. Also let X (Y) and a (o>) be the mean and x monthly mean thermal atmospheric forcing, but with standard deviations of the reference (estimate) values, wind forcing that includes the 6-h sub-monthly varia- respectively. The statistical relationships (e.g., Murphy, bility because of mixed layer sensitivity to high 1995) between NOAA and HYCOM SST time series at a given grid point are expressed as follows: frequency forcing (e.g., Kara ct al.. 2005a). Performing a 1-year global HYCOM simulation using ME= Y -X. (6) any of the atmospheric forcing products required ~ I 5 wall-clock hours on 64 HP/Compaq SC45 processors. t 1/2 Thus, the 0.72° global HYCOM provides inexpensive, x; (7) non-cddy-resolving ocean model simulations. RMS Monthly mean SST fields obtained from the model simulations arc evaluated through extensive model-data comparisons using various statistical metrics (Section R = -J2 {Xi-X)(Yi-Y)/(a a ), (8) 4.1). Within the quantitative framework, statistical error x r and skill analyses are then presented for quality assessment of the model in simulating climatological 2 2 (9) monthly mean SST using bulk parametcrizations and the SS = R -[R-(<T /<rx)} - [(Y X)/a \ Y x three atmospheric forcing products (Section 4.2). ME (i.e., annual bias) is the mean error between the 4.1. Statistical metrics mean HYCOM and NOAA SST values. RMS (root- Monthly mean HYCOM SST climatologies are mean-square) SST difference is an absolute measure of constructed from SSTs obtained from the last 4 (out of the distance between the two time series, and the R 8) model years. For example, the climatological mean value is a measure of the degree of linear association January SST is formed using mean January SSTs from between the time series. The non-dimensional SS given in Eq. (9) is the fraction of variance explained by model years 5-8. The same process is repeated for the other months. Because all simulations performed with HYCOM minus two dimensionless biases (conditional the 0.72° model use climatological mean forcing from bias, # |, and unconditional bias, # , .„ i) which arc com u 1( m NWP products and there are no mid-latitude mesoscale not taken into account in the R formulation (8) as explained in Murphy (1988), in detail. B , described eddies, the 4-ycar averaging period is sufficient. uncand For model-data comparisons the NOAA SST clima- as systematic bias, is a measure of the difference between tology (Reynolds ct al., 2002) is taken as a reference the means of NOAA and HYCOM SST time series. i? ] mni (truth) because its resolution (l°x 1°) is close to that of is a measure of the relative amplitude of the variability in IIYCOM (0.72° * 0.72° cos(lat)). The NOAA SST fields the NOAA and HYCOM time scries or simply a bias due