Boemre 2008 Extended Hindcast Calculation of Gulf of Mexico Circulation: Model Development, Comparisons with Observations



Download 265.82 Kb.
Page2/4
Date18.10.2016
Size265.82 Kb.
#894
1   2   3   4

1. INTRODUCTION
The Gulf of Mexico is a semi-enclosed sea in the southwestern tropical/subtropical Atlantic Ocean (fig.1.1).  The Gulf has two narrow openings.  The Yucatan Channel is V-shape between the Yucatan Peninsula (west) and Cuba (east); it is approximately 200 km wide at the surface and 25 km at its sill depth of approximately 2040m (Sheinbaum et al., 2002), and connects the Gulf to Caribbean Sea.  The Straits of Florida is between Cuba (south) and Florida (north); it is approximately 100 km wide at the surface, has a shallower sill depth of about 800m, and connects the Gulf to the Atlantic Ocean.
fig.appendix.1_20100122edav3_lyo.tif

Figure 1.1 The northwest Atlantic Ocean model (  6~12 km within the Gulf of Mexico) (NWAOM; the whole region shown) and the fine-grid resolution (  3~7 km) Gulf of Mexico and northwest Caribbean Sea region (west of the dashed line). Contours show isobaths in meters, and silhouettes at 55oW show steady transports specified for the Gulf Stream and returned flows.


The Yucatan Current is a narrow and swift western boundary current (WBC) along the eastern coast of Yucatan Peninsula in the northwestern Caribbean Sea (the Cayman Sea). Its mean speed is about 1 m s-1 but at any time this can exceed 2 m s-1 (Cetina et al. 2006). The Yucatan Current supplies approximately 23 to 27 Sv (1 Sv = 106 m3 s-1; Schmitz et al. 1996; Candela et al. 2002) of warm and relatively more saline waters from the Caribbean Sea into the Gulf of Mexico (Nowlin, 1972). Temperatures T (potential) in the near-surface 100~200 m are generally above 26oC with salinities S of about 36.4 psu. At T  6.3oC in the subsurface depths z  1000 m the salinity reaches a minimum < 34.9 psu that is a characteristic of the remnant of Subantarctic Intermediate Water (Wϋst, 1964). In the upper 500~1000 m, the Yucatan Current intrudes northward into the Gulf along the Campeche Bank’s eastern shelfbreak, and loops clockwise into the Straits of Florida. This strong clockwise-turning current in the eastern Gulf of Mexico is called the Loop Current. The Loop Current is highly inertial and, because of its generally large size of  300 km or larger in its south-north and west-east extents, is significantly affected by the variation of the Coriolis parameter (f), i.e. by the beta () effects (Hurlburt and Thompson, 1980; Nof, 2005). Large warm-core rings (Loop Current eddies, 200~350 km wide and 500~1000 m deep) episodically separate from the Loop Current at time intervals (“periods”) that range from a couple of months to almost two years (Sturges and Leben, 2000). An idealized but useful idea of why eddies are shed is given by Nof (2005). The Loop Current grows with mass influx from the Yucatan Current. Up to a certain large size, the westward Rossby wave speed (which is  R2, where R is the Rossby radius based on the matured deep Loop) overcomes the growth rate and the Loop sheds an eddy. Because of their dominance, the Loop Current and Loop Current eddies affect just about every aspect of oceanography of the Gulf of Mexico (Oey et al. 2005).
The threat of an oil spill contacting land in the Gulf of Mexico is of great concern. The 2010 BP oil spill is a vivid example [Chang et al. 2011]. To estimate the potential for such a spill coming into contact with resources in the region, a robust knowledge of the circulation that accounts for the dominant physical processes in the Gulf is needed. Moreover, extensive surface and subsurface current data such as those produced over years of observational and modeling efforts sponsored by BOEMRE are crucial in oil spill prediction and projection [see, e.g. Chang et al. 2011]. Oey et al. [2005a] provide a comprehensive review of these processes which beyond the shelfbreak (water depths deeper than approximately 500 m) are generally dominated by powerful features with scales ranging from a few tens (in small eddies and fronts) to hundreds (in warm rings) of km’s, and also by the Loop Current in the eastern Gulf [Hamilton and Lugo-Fernandez, 2001; Hamilton, 2007; Oey, 2008; Hamilton, 2009]. The near-surface layers of the Gulf, as well as the shelves, are also affected by atmospheric forcing such as the winter storms and hurricanes during the summer and fall [Oey et al. 2006; 2007]. Over the shelves, buoyancy forcing from rivers is also important [Oey, 1995]. Strong winds produce not only intense currents but also large waves [Wang and Oey, 2008; Mellor et al. 2008], and there is also strong interaction between mesoscale eddy-features and storm-induced currents that can produce intense subsurface currents [as deep as 500m below the surface: Oey et al. 2008]. To produce a comprehensive current dataset that can mimic as close as possible the state of the real ocean, a time-dependent three-dimensional circulation model assimilated (or “injected”) with observational data is deemed necessary [e.g. Oey et al. 2005b; Lin et al. 2007; Yin and Oey, 2007]. Such an effort is attempted here, for the period 2000 to 2007.
3.1 Background of Princeton Modeling Efforts
Princeton University has previously conducted numerical modeling of the Gulf of Mexico for BOEMRE [Oey, 2004]. The model domain encompasses the entire US east coast, the Gulf of Mexico in particular (Fig.1.1; the northwest-Atlantic Ocean model or NWAOM domain). The simulation covers the year 1993-1999, and uses as the forcing (1) six-hourly winds from ECMWF analysis, (2) monthly surface mass fluxes (corrected by satellite SST, see below), (3) daily discharges from 34 US rivers, (4) climatological density fields as initial conditions and also to prevent long-term drift, and (5) climatological density and transport at the eastern model boundary in the Atlantic Ocean (Fig.1.1). To correctly simulate the positions of major (and very important) oceanic features such as the Loop Current and rings in the Gulf, an efficient data-assimilation scheme (based on optimal interpolation, OI) is used whereby satellite altimetry and sea-surface temperature data are combined with the model dynamics to produce a hindcast field (sea surface height, currents, temperature, salinity and other diagnostic variables such as the turbulence eddy viscosities and diffusivities). Data assimilation is the process of combining a physical model with observational data to provide a state analysis of the system which is better than could be obtained using just the data or physical model alone. Figure 1.2 shows an example of such data-assimilated current (vectors) and sea-surface height (red is high and blue is low) fields. There is visual agreement between the model and satellite (black contours) positions of Loop Current and rings. In fact, when compared to independent in situ (at the “” location shown in Fig.1.2) and ship data, the model fields are more accurate than satellite products especially at smaller scales. Mean speed and direction (absolute) errors are 0.4~5 cm s-1 (1~10% errors) and 10o~20o respectively, and amplitude and phase of model-to-observed complex (velocity) correlation are 0.76~0.86 and 0.3o~7o respectively [Lin et al. 2007]. We have also conducted other extensive checks with our model, both for its accuracy when compared against observations, and also for its ability to correctly simulate processes (dynamics; see publication list in http://aos.princeton.edu/WWWPUBLIC/PROFS/).

fig02_velmdrf_elmsat_hcast2me_pseudo_ssh_nx=5

Figure 1.2. Daily-averaged analysis sea-surface height (SSH; color) and surface currents shown on the first of each month from Sep/1999 through Aug/2000 in the Gulf of Mexico. Thick dark contour is the zero-value of satellite SSH. [From Lin et al. 2007].


Subsequent studies supported by BOEMRE further improve the above model since it and its results for the period 1993-1999 were first delivered to BOEMRE for Contract#1435-01-00-CT-31076 (COTR: Dr. Walter Johnson). We have incorporated a nesting scheme [based on Oey and Chen, 1992] that doubles the model’s horizontal grid resolution “” in the Gulf (to become  3~5km), and at the same time the nest is forced by the parent grid and retains all functionality and assimilative capability of the parent grid [Fig.1.1; Oey and Zhang, 2004; Oey, 2007]. We have also changed our wind field to a higher-resolution wind dataset (1/4o×1/4o instead of 1o×1o) that also has satellite data wind blended into it, and moreover an algorithm has been developed to correct the wind input for intense wind field due to tropical storms using high-resolution and accurate analysis from the NOAA’s Hurricane Research Division (HRD). The procedure is such that different storms can follow one another or storms can co-exist during some overlapping period; see http://aos.princeton.edu/WWWPUBLIC/PROFS/ANIMATION/hwind/hwind_movie_2005_gfs.html for an example animation. The model also now has improved mixed layer physics that account for the effects of (wind) wave-breaking. All these features are incorporated in the model hindcast conducted for this project. We also incorporate two additional features: a numerical scheme that further reduces the truncation error in regions of steep topography, and a data assimilative scheme based on the Ensemble Kalman Filter methodology. The various methods will be evaluated in the following sections.
Observational Data
Additional observations have become available since the above-mentioned “1993-1999” hindcast study was completed in Oey [2004]. They are: the Exploratory Study of Deepwater Currents in the Gulf of Mexico, the Survey of Deepwater Currents in the Western Gulf of Mexico, and Survey of Deepwater Currents in the Northeastern Gulf of Mexico. There are also data from the industry, e.g. the “MMS ADCP”1 data on the NDBC website (since May/2005). There are also data from the LSU moorings, the “Explorer” data (Yucatan Channel and Florida Straits), the Mexican “DeepStar” data across the Yucatan Channel, SAIC’s Sigsbee moorings, and the PALACE float data (at 900m). We plan to use some of these observations for independent skill-assessment, or for direct assimilation into the model. The accuracy of hindcast results will be evaluated in details.

2. THE MODEL
The terrain-following (i.e. sigma) coordinate and time-dependent numerical model for this study is based on the Princeton Ocean Model (Mellor, 2002). The Mellor and Yamada’s (1982) turbulence closure scheme modified by Craig and Banner (1994) to effect wave-enhanced turbulence near the surface is used (Mellor and Blumberg, 2004). A fourth-order scheme is used to evaluate the pressure-gradient terms (Appendix 1; see Berntsen and Oey, 2010) and, in combination with high resolution and subtraction of the mean -profile, guarantees small pressure-gradient errors of O(mm/s) (c.f. Oey et al. 2003). The Smagorinsky’s (1963) shear and grid-dependent horizontal viscosity is used with coefficient = 0.1, and the corresponding diffusivity is set 5 times smaller (c.f. Mellor et al., 1994). Two domains are used on orthogonal curvilinear grids (fig.1.1). The outer grid is the regional northwestern Atlantic Ocean model (98W-55W and 6N-50N; hereinafter referred to as NWAOM) reported in our previous studies (e.g. Oey et al. 2003). The NWAOM has 25 vertical sigma levels and horizontal grid sizes   6~12 km in the Gulf of Mexico. The World Ocean Atlas data (“Climatological” data) from NODC (http://www.nodc. noaa.gov/OC5/WOA05/pr_woa05.html) was used for initial condition as well as boundary condition along the eastern open boundary at 55oW. Across 55W, a steady transport combined with radiation using also the geostrophically-balanced surface elevation g (Oey and Chen, 1992a) specifies the Gulf Stream exiting near the Grand Banks south of Newfoundland with a magnitude of 93 Sv following Schmitz (personal correspondence, see also Schmitz, 1996; Hendry 1982; Hogg 1992; Hogg and Johns 1995). This is balanced by transports specified as broad return flows south (the “Worthington Gyre” - Worthington, 1976) and north (the “North Recirculation gyre” - Hogg et al. 1986) of the jet. The vertical structures of the currents (i.e. after a depth-averaged value is removed) are specified using radiation conditions. The velocity component tangential to the boundary, as well as turbulence kinetic energy and length scale, are specified using one-sided advection scheme at outflow grids and are set zero at inflow. The (potential) temperature (T) and salinity (S) are similarly advected during outflow, but are specified using climatological values at inflow grids. Radiation is used for the surface elevation , but since POM uses a staggered C-grid, and because transports are specified, the boundary  plays only a minor role and a zero-gradient condition on it works well also. At the sea surface, heat and salt fluxes are specified using a combination of monthly climatological values from NODC [http://www.nodc.noaa.gov/OC5/WOA05/ pr_woa05.html] and relaxation to climatological sea-surface temperatures and salinities (see Oey and Chen, 1992a). The model is also forced by 6-hourly wind (described later), as well as by daily discharges from 34 rivers along the northern Gulf coast according to the method given in Oey [1995; 1996]. To prevent temperature and salinity drift in deep layers in long-term integration, the T and S for z < 1000 m are (weakly) restored to monthly-mean climatological values with a time scale of 600 days. More details are in Oey et al. (2003) and Oey (2004). The NWAOM is continuously run from 1992 through present with data assimilation (next section). The 5-day averaged results are used to obtain boundary values for the second, nested domain, described next.
The second domain (98W-79W and 8N-30N; Fig.1.1) contains the Gulf of Mexico and a portion of the northwest Caribbean Sea (i.e. the Cayman Sea). It has a doubled-resolution grid (  3~7 km, same 25 sigma levels). It too has only one open boundary in the east (at 79W) and the boundary specifications are similar to those described above for the NWAOM, except that for T and S a flow-relaxation scheme with 20 cells for the relaxation zone (Oey and Chen, 1992b) is used near the Caribbean (inflow) portion of the boundary. Also, the NWAOM values (time-dependent 5-day averaged transports, T and S) are used instead of climatological data. This one-way nesting follows Oey and Zhang (2004). The high-resolution nested domain is then run from 1992-2007, also with data assimilation, with 6-hourly wind, with daily discharges, with monthly surface heat and evaporative fluxes and with weak restoring for layers deeper than z = 1000 m.
2.1 Wind Forcing
Wind is an important part of the input forcing to the model. This is so in particular for oil spill trajectory assessments near the surface. Therefore, care must be taken to ensure that the data used is the best available. For the present project, we use the Qscat/NCEP blended wind which is available every six hours on a 0.5o×0.5o grid. This wind is read and temporally and spatially interpolated “on the fly” (as the model is running) onto the model domain. The wind velocity is then converted to wind stress. For a limited number of runs for research and testing purposes, for example in the oil spill trajectory calculations, the “Cross-Calibrated Multi-Platform (CCMP) wind products were also used [Chang et al. 2011]. The CCMP was not available until after the calculations for the present project ended and was therefore not used in the calculations that were delivered to BOEMRE as part of the required deliverables. However, for both products, we conducted extensive comparisons, both amongst themselves as well as between them and measured winds at NDBC buoys. The CCMP is generally more accurate, and while there are some instances when the two datasets differ, both achieve high vector correlations when compared with the buoy winds, with values generally larger than about 0.85 and small deviation (error) angles of only a few degrees.

fig2.1_final_report.png

Figure 2.1. Bi-monthly variance ellipses and mean (vectors originating from ellipses’ centers) for wind stress derived from the Qscat/NCEP blended product, computed for 1988-2008. Vector scale (0.2 N m-2) is shown for the mean wind stress, and color background indicates the standard deviation, also in N m-2. Very similar plots are obtained using the CCMP wind.


To calculate wind stresses, we use a bulk formula with a high wind-speed limited drag coefficient that fits data for low-to-moderate winds (Large and Pond, 1981) and data for high wind speeds (Powell et al. 2003), as given in Oey et al. (2006):
Cd 103 = 1.2, |ua|  11 m s-1;

= 0.49 + 0.065 |ua|, 11 < |ua|  19 m s-1;



= 1.364 + 0.0234 |ua|  0.00023158 |ua| 2, 19 < |ua|  100 m s-1 (1)
where |ua| is the wind speed.2 According to this formula, Cd is constant at low winds, is linearly increasing for moderate winds, reaches a broad maximum for hurricane-force winds, |ua|  30~50 m s-1, and then decreases slightly for extreme winds. Donelan et al. (2004) suggest that the Cd-leveling at high wind may be caused by flow separation from steep waves. Moon et al. (2004) found that Cd decreases for younger waves that predominate in hurricane-forced wave fields. Bye and Jenkins (2006) attribute the broad Cd-maximum to the effect of spray, which flattens the sea surface by transferring energy to longer wavelengths.
Figure 2.1 shows the bi-monthly mean and variance ellipses computed from the CCMP product for 1988-2008 (the Qblend product is very similar; see Appendix 3). It shows the seasonal transition of the mean wind and its variability in the model region. The strongest mean wind and variance are during late fall through early spring (figs.2.1f,a,b) concentrated along the storm track coincident with the Gulf Stream (e.g. Minobe et al. 2008). In the Gulf of Mexico, the wind is weaker. The mean wind is easterly and strongest in winter with peak mean wind stress of approximately 0.7 N m-2 (c.f. Gutierrez and Winant, 1996). In summer the wind becomes southerly and southeasterly. Variances are larger than means. Figure 2.1e also shows large variance in the northeastern and central Gulf in September and October as a result of tropical cyclone activity during summer and early fall. Wind animations from 1999 through 2008 are given as *gif files in ftp://aden.princeton.edu/pub/lyo/mms/walt/mab/qnblend_windstrs/.
The Qblend wind is still inadequate to account for intense hurricane winds. This is due both to the lack of spatial and temporal resolutions. So we combine the Qblend wind with the Hurricane Research Division (HRD; http://www.aoml.noaa.gov/hrd/) wind. The method was previously described in Yin and Oey (2007). The HRD data is given in a 1000 km  1000 km (dimensions are approximate) moving “box” centered about the hurricane’s track. Storm centers are first linearly interpolated to hourly locations. Consecutive HRD maps are then overlapped at the hourly locations and linearly interpolated. The hourly maps are then merged with Qblend wind using a weight that retains the HRD data within a circle of radius = 0.8side of the box (i.e.,  400 km), and that smoothly (using a tanh function) merges the HRD and NCEP winds beyond that radius. Figure 2.2 shows an example of the merged winds for Hurricane Katrina (August 24-29, 2005).

fig04_hrd_gfs_wind_plot_gfs_aiedited

Figure 2.2. Twelve hourly plots of HRD/NCEP winds showing the path of Hurricane Katrina from (A) August 24 at 12:00 GMT through (K) 29 at 12:00 GMT, 2005. The last panel (L) is for August 29 at 19:00 GMT. Dots indicate daily locations of the storm’s eye.
2.2 Baroclinic Pressure-Gradient Error at Steep Topography
Sigma-coordinate ocean models such as the one we are using are attractive because of their abilities to resolve bottom and surface boundary layers. However, these models can have large internal pressure gradient (IPG) errors. Oey et al. (2003) have shown that the IPG errors for our model in the region shown in fig.1.1 are small. However, model bias should be carefully examined when interpreting physical processes in regions with steep topographic gradients. In Berntsen and Oey (2010; which is included in Appendix 1), two classes of methods for reducing the IPG errors are assessed. The first is based on the integral approach used in the Princeton Ocean Model (POM). The second is suggested by Shchepetkin and McWilliams (2003) based on Green’s theorem; thus, area integrals of the pressure forces are transformed into line integrals. Numerical tests on the seamount problem, as well as on a northwestern Atlantic grid using both classes of methods, are presented. For each class, second-, fourth-, and sixth-order approximations are tested. Results produced with a fourth-order compact method and with cubic spline methods are also given. The results show that the methods based on the POM approach in general give smaller errors than the corresponding methods given in Shchepetkin and McWilliams (2003). The POM approach also is more robust when noise is added to the topography. In particular, the IPG errors may be substantially reduced by using the computationally simple fourth-order method from McCalpin (1994). Stream function error maps for a grid resolution that is 2 times coarser (hence the error is 4 times larger) than that used in the outer NWAOM grid of fig.1.1 (Oey et al. 2003; also Oey, 2004) are given in fig.2.3. The forth-order methods (McCalpin4th in fig.2.3b and compact differencing CompactD in fig.2.3c) both give small errors O(0.1~0.3 Sv) in the Gulf of Mexico. For the NWAOM grid, the corresponding error is O(0.25~0.75 Sv), and for the finest nested grid it is O(0.06~0.19 Sv). These are much smaller than the O(10~30 Sv) in the Loop Current and rings, and O(1~3Sv) which may exist in the deepest regions of the Gulf.

Figure 2.3 Stream functions (+ for cyclonic circulation) after 180 days for a the pom2nd method, b the pom4thMcC method, c the pom4thCmp method, d the sm4thCubicH method, e the sm4thpoly method, and f the pom6th method. For clarity, small values for water depths less than 200m are omitted.


3. DATA ASSIMILATION
Loop Current and associated rings dominate the upper-layer (z  0 to 1000m) circulation in the Gulf of Mexico [e.g. Oey et al 2005]. A recent collection of papers in a book edited by Sturges and Lugo-Fernandez [2005] give useful information on the physical oceanography of the Gulf of Mexico. For modelers and those interested in ocean predictions involving meso-scale features, the Loop Current and rings are “ideal.” Some of the reasons are because the Loop Current and rings (1) are energetic (speeds 1~2 m/s), hence have large signal-to-noise ratio; (2) have fairly large scales (> 200 km), hence are more easily resolvable by most modern models with grid sizes of approximately 10 km and less, and also by satellites; (3) are entirely contained within the Gulf of Mexico  a basin of large-enough scales (~3000 km) that dynamics can be affected by planetary beta effects; (4) can be significantly affected by their interactions with slopes and shelves; and (5) offer excellent case studies of interaction between energetic ocean currents (eddies) and tropical cyclones that frequently visit the Gulf during the hurricane season. On the other hand, the same reasons listed above contribute also to complex dynamics that make it challenging to simulate and predict the Gulf’s circulation; in particular, the behaviors of the Loop Current and rings.
Assimilating observations into models is one way to constraint the solution to derive an initial model state that is closer to the observed, presumably the “truth” state. However, even with a good initial state, the strong nonlinearity almost certainly guarantees that the simulated Loop Current and rings will, after some time, depart from the (future) observed state. Using the breeding method [e.g. Kalnay, 2003], Yin and Oey [2007] introduce a random field of initial perturbations into their primitive-equation circulation model (the Princeton Ocean Model; Mellor, 2004) of the Gulf of Mexico, and show that a fully-evolved bred vector consisting of cyclones and anticyclones develop in approximately 6~8 weeks. These are also the time scales in which their ensemble model forecasts tend to depart significantly from satellite altimeter observations of the positions of the Loop Current and rings. Similar conclusions are also obtained by Oey et al [2005]. The authors focus on predicting the detailed frontal evolutions of the Loop Current and rings, and show that one large source of error came from the inexact knowledge of the positions of these fronts. It is clear that good data-assimilative analyses are necessary for good forecasts.
The assimilation schemes include the Mellor and Ezer’s [1991] statistical correlation and optimal interpolation method as well as the Ensemble Kalman OI [Evensen, 2006] methods, and tests with various model grid resolutions both in the horizontal and vertical are also conducted. In this work, we assimilate satellite SSHA (o, from AVISO, www.aviso.oceanobs.com; Ducet et al. 2000) and SST (from the United States GODAE, www.usgodae.org) to derive the model hindcast analysis.
3.1 Statistical Optimal Interpolation Scheme
Satellite data are assimilated into the model following the methodology given in Mellor and Ezer (1991) and Ezer and Mellor (1994). In this method, SSHA is projected into the subsurface temperature field using pre-computed correlation factors derived from a long-time ( 10 years) prognostic integration that has yielded a statistical equilibrium eddy field. Thus the resulting temperature anomaly (T) is (<.> is time-averaging, and T is the potential temperature):
T(x,y,z ,t) = FT(x,y,z) o(x,y,t), (3.1)
where the correlation factor is ( = model SSHA)
FT = <T >/<2>, (3.2a)
and the corresponding correlation coefficient is
CT = <T >/(<T2><2>)1/2. (3.2b)
Ezer and Mellor (1994) assimilate along-track o data assuming a linear-saturation error growth model for the first-guess error. Our experience has been that if AVISO o maps are assimilated the following simplified formula (see Wang et al. 2003) suffices:
Ta = T + (2 RACT2/(1 + 2 RACT2  CT2)) (TO  T) (3.3)
where T is the model (first-guess) temperature, Ta denotes the analysis temperature, RA is the ratio of the assimilated time step tA to the de-correlation time scale tE of the model eddy field, and TO is the ‘observed’ temperature inferred from (3.1),
TO = + FT o. (3.4)
Instead of using the model mean for in (3.4), our past experience has been that setting = TC, the observed temperature climatology, helps to control long-term (~ 10 years) drift in the model. For the present application, the differences are small. Formula (3.3) assumes that the AVISO map errors are small compared to the model errors, and that tA << tE. We follow Ezer and Mellor (1994) and set tA = 1 day. The tE is estimated from the above-mentioned 10-year prognostic model run and is  30 days in regions of the Gulf of Mexico dominated by the Loop Current and rings. This may be compared with the value of 20 days used by Ezer and Mellor’s (1994) for the Gulf Stream which therefore appears to have shorter meander and eddy evolution time scales. The tE is also proportional to the time scale of the model error growth, and the 30-day value is consistent with Oey et al’s (2005b) findings of predictability time scales of about one month for the Loop Current and its associated rings. As pointed out by Ezer and Mellor (1994), the assimilation (3.3) is such that Ta  TO in regions where the correlation is high (CT2  1), but Ta  T where the correlation is low. A similar assimilation of SST is also carried out after (3.3) with CT and FT replaced by the corresponding functions that use (SST) in place of  in (3.2). The SSHA and SST assimilations complement each other: SSHA assimilation is most effective over deep waters (for isobath > 500 m) while SST assimilation influences waters on shallow shelves. For more details see Wang et al. (2003), Fan et al. (2004), and Oey et al. (2005b).
3.2 The Ensemble Kalman Filter Scheme:
A sequential data assimilation method is used:

, (5a)

, (5b)

where is forecast state, is observation, is analysis state, is observational operator, is gain matrix, is forecast error covariance, is observation error covariance, denotes transpose, is matrix inverse. We assume that is known. The ensemble optimal interpolation (EnOI) method is used for assimilation; the error covariance is not fixed however (below). Instead of equation (5b), we use the following formula to calculate :



, (6)

where is a parameter that controls the weight between and .


3.2.1 Estimating the Forecast Error Covariance:

The forecast error covariance (Pf) is estimated by an ensemble method. In the full Ensemble Kalman Filter (EnKF) method, the Pf is time-dependent calculated at every assimilation time-step ta (= 1 day, say). One way is to compute at every ta the Pf from the solutions of a set of perturbed experiments. In EnOI, the Pf is fixed; one way is to take the temporal mean of error covariances from a multi-year integration, and then use that mean as the time-independent Pf in (5) and (6) for the analyses covering the same multi-year period. The EnKF is very cpu-intensive; the EnOI is simple but it does not take into account the “error of the day.” In this work, we have designed a slightly different algorithm that partially takes into account the evolution of the background field. Instead of computing the Pf at every ta as done in EnKF, the ensemble of perturbed experiments are conducted every tE, where tE corresponds to the time scale in which the eddy field has evolved “significantly.” Yin and Oey’s [2007] bred-vector analysis shows that random perturbations would evolve into well-organized eddy field in approximately 8 weeks. This time scale agrees with past experiences of modeling the Loop Current and rings [e.g. Oey et al. 2003]; we take, then, tE = 60 days.


The ensemble is comprised of, then, a set of model solutions corresponding to the analysis time interval ta. The initial states of different ensemble members are given small random perturbations [see Yin and Oey, 2007 for details], and these ensemble calculations are done every tE = 60 days. There are ensemble solutions denoted as (), where is a column vector including all model variables (, etc). Then the forecast error covariance is estimated by

, (7)

where matrix , vector .



contains the correlation between any two model variables, as well as between two model grids of one model variable.
Theoretically, the estimated will approach the true as approaches infinity. However, this is not practical, so there is likely some spurious covariance in at large spatial distances. The following piecewise continuous function (Gaspari 1999) is adopted to filter out this kind of spurious covariance.

, (8)

where is a given positive constant, is the input distance.


In summary, then, our assimilation method is a little more sophisticated than the standard EnOI method. It adopts the EnKF idea in that the changes with time, though discretely at 60-day intervals. The is calculated by ensemble experiments with 30 ensemble members each of which is initially given random perturbations. We will refer to our scheme as the QEnKF – Quasi-EnKF.
3.2.2 Assimilating Altimetry SSHA & Moored Temperature Data
We first describe how a two-dimensional (xy) dataset such as the altimetry SSHA is assimilated. To assimilate at the model grid , the following algorithm is used:


  1. Given a small (xy) region containing the model grid in a circle (say of radius < 1.0 degree; Figure 3.1a), find all the SSHA observational grids inside this region;




  1. For each , find the smallest quadrilateral that encloses it, formed by four model grids , , , ; one of these model grids is - the point to be assimilated;



  1. Use all , , , and to define the observational operator ;



  1. Compute , the covariance of SSHA on observational grids . Compute = the covariance between temperature model grid and the SSHA on observational grids .



  1. Compute gain matrix using and ;



  1. Compute analysis temperature at .

To assimilate the moored temperature observations, the following algorithm is used (see Figure 3.1b):
(a) Given a small (3-D) region within the model domain, find all the moored temperature observations inside this region;
(b) For each , find the smallest quadrilateral with four model grids in horizontal direction: , , , . In addition, find two model vertical level and such that . Then we obtain a 3-D grid box with eight model grids:

, , , and , , , .

(c) Use all above , , , , , , , and


to derive the observational operator H.

(d) Compute , the covariance of temperature observational grids . Compute also , the covariance between temperature model grid and the temperature observational grids .


(e) Compute gain matrix using above and .
(f) Compute analysis temperature for .
fig02_sun_oey_da_lce_rings_deep_v01.tif

Figure 3.1a Illustrating how observational points are defined within the neighborhood of a model grid.


fig03_sun_oey_da_lce_rings_deep_v01.tif

Figure 3.1b Illustrating how temperature (or any variable) is assimilated at subsurface.





Download 265.82 Kb.

Share with your friends:
1   2   3   4



The database is protected by copyright ©ininet.org 2024
send message
    Main page
gulf stream
scientific community
northern gulf
main text
eastern gulf
loop current
region north