The standards for skill assessment of operational marine

Chinese Journal of Oceanology and Limnology Vol. 25 No. 1, P. 27-35, 2007 DOI: 10.1007/s00343-007-0027-7

The standards for skill assessment of operational marine forecast system*

ZHANG Aijun (张爱军) , ,**, FAN Wenjing (范文静) , JI Fengying (纪风颖)

( National Marine Data and Information Service, Tianjin 300171, China)

( Coast Survey Development Laboratory, NOS/NOAA, 1315 East-West Highway, Silver Spring, MD 20910, USA)

Received Aug. 8, 2005; revision accepted June 16, 2006

Abstract To support navigational and environmental applications in coastal waters, marine opera-tional forecast models must be developed and implemented. A forecast model must guarantee that it is scientifically sound and practically robust for performance and must meet or excel all target frequencies or durations before being released to the public. This paper discusses the standard policies and procedures for evaluation of operational marine forecast models. The primary variables to be evaluated are water lev-els, currents and water density (water temperature and salinity). Key words: skill assessment; hydrodynamic models; forecast

1 INTRODUCTION

An operational ocean forecast system is de-signed to provide information regarding both pre-sent (nowcast) and future (forecast) ocean condi-tions at many locations such as estuary, bay, or coastal area. The primary variables of the system are water level, current speed, salinity, and temperature. Some key products derived from an operational ocean weathercast system are: (1) nowcasts that provides the most useful and accurate description of present state of the sea including living resources; (2) forecasts that predicts future condition of the sea as far ahead as possible; (3) hindcasts is to assemble long term data sets for providing data descriptive of past weather states, and time series showing trends and changes. The processed data as final products must be delivered rapidly to customers in industry, government, and other organizations including regulatory authorities to help protecting marine en-vironment, improving safety operation in marine transport and construction, warning public risk in emergent situations against floods, tsunami, tropical storm, and so on.

All forecast systems must be skill-assessed be-fore they are put into practice as operational forecast systems to ensure that all models have been devel-

oped and implemented in a scientifically sound and operationally robust manner. This paper discusses the specified procedures for evaluating an opera-tional ocean forecast model. Skill assessment is an objective measurement of how well the model nowcast or forecast performs when compared to observations. The approach here is to test the per-formance if it can (1) simulate astronomical tidal variability; (2) simulate total (tide and non-tidal ef-fects) variability in both the model development stage and the real-time operation, and (3) yield a more accurate forecast than the tide tables and/or persistence. A standard skill assessment scores automatically using data files containing observed, nowcast, and forecast variables. These data are processed for skill assessment are displayed in ta-bles into which model evaluation reports incorpo-rated. The procedures include harmonic analysis, gap filling, filtering (or singular value decomposi-tion), and extracting water level and current ex-trema.

The data time series required by skill assess-ment are defined in Section 2. The techniques of data analysis for processing and analyzing data se-ries are introduced in Section 3. The standard pro- * Supported by the National Natural Science Foundation of China (No. 40376010).

** Corresponding author:

28CHIN. J. OCEANOL. LIMNOL., 25(1), 2007 Vol.25

cedure of skill assessment is discussed in Section 4

later.

2 DATA REQUIREMENTS

Skill assessment for each location requires 3 basic types of time series data: observed, tidally predicted (for tidal regions), and model simulated. A uniform time interval of 6 min is required for each series, but 1 h interval is suitable for water levels. The ideal duration of each time series is 365 days in order to capture all expected seasonal conditions. However, it is sometimes difficult to get such a long time series. Therefore, the minima of 6 months for water level and 29 days for current are recom-mended. All model and observational data units must conform to the International Standard for units and time reference (UTC).

All observational data should be quality con-trolled and processed to final units (e.g., meter or m/s). It is expected that there will be occasional gaps that can be filled by some simple methods Tidally predicted data are based on 37 standard con-stituents obtained either from the accepted harmonic constants or derived from observational time series by harmonic analysis. The standard prediction method uses harmonic constants, lunar node factors, and equilibrium arguments.

Model simulated data are generated by running the model under one of 4 scenarios: astronomical tide only, hindcast, semi-operational nowcast, and semi-operational forecast. The scenarios are de-scribed below.

2.1 Model-running scenarios

There are 4 scenarios under which the model is run to produce data for the skill assessment, and they are discussed in the order of occurrence during model development.

2.1.1 Astronomical tide simulation only

For regions where significant tidal variations exist, the model is run in astronomical-tide-only scenario as tidal variations may account for a sig-nificant part of error. In this scenario, the model is forced with harmonically predicted astronomical tides only for ocean boundary water levels with no surface forcing (wind, pressure, etc.). The tempera-ture and salinity are set a constant with zero river flow input. The model time series can be compared

with tidal predictions, and harmonically analyzed to produce constituent amplitudes and phases for comparison with accepted values. 2.1.2 Hindcast

The model forcing is based on historical, best available gap-filled observational data for open boundary water levels, surface winds, temperature, salinity, and river flows. The model time series can be compared with the available observations. 2.1.3 Semi-operational nowcast

The model forcing is based on real time ob-served values. The real-time observation may be incomplete with gaps. The operational model will be restarted more often (for instance, 4 times a day). The ability of the model to work correctly in the re-starting mode will be tested. This run tests the abil-ity of the model under real working conditions. 2.1.4 Semi-operational forecast

The model forcing is based on recent forecast values from other models, even though some data could be missing. Initial conditions are generated from observed data or the output from a nowcast. This run tests the ability of the model under working conditions.

2.2 Time series variables

The following time series are required for skill assessment computations. The definitions are sum-marized in Table 1.

There are 3 groups of data sets.

Group 1: the data can be either a time series of values (such as observations at a location) or a se-ries of values from concatenated segments (such as a set of 24 h nowcasts or forecasts starting at one time in the day). For currents, speed and direction will be needed for the time series; the direction error is computed only for current speeds greater than 0.26 m/s.

Group 2: a set of values are created from Group 1 series by selecting a sub-set of values such as the time and amplitude of high water or the time of the start and end of slack water (defined as hav-ing a current speed less than 0.26 m/s).

Group 3: the data consist of the values of the forecast variable that are valid at a fixed interval into the forecast (e.g., 0, 6, 12 h, etc). The compari-son series is then the observed variable at the time

No.1 ZHANG et al.: The standards for skill assessment of operational marine forecast system

Table 1 Data series groups and the variables in each group*

Group Variable description

Group 1 (Time series)

Water level

Current speed Salinity

Group 2 (Values at a tidal stage)

Group 3 (Values from a forecast)

*

29

Symbol H, h U, u S, s AHW, ahw ALW, ahw THW, thw TLW, lw AFC, afc AEC, aec TFC, tfc TEC, tec DFC, dfc DEC, dec TSF, tsf TEF, tef TSE, tse TEE, tee Hnn, hnn Unn, unn Dnn, dnn Snn, snn Tnn, tnn

Current direction D, d Water temperature T, t

Amplitude of high water Amplitude of low water Time of high water Time of low water

Amplitude of maximum flood current Amplitude of maximum ebb current Time of maximum flood current Time of maximum ebb current Direction of current at maximum flood Direction of current at maximum ebb Time of start of current slack before flood Time of end of current slack before flood Time of start of current slack before ebb Time of end of current slack before ebb Water level at forecast projection time of nn hrs Current speed at forecast projection time of nn hrs Current direction at forecast projection time of nn hrs Salinity at forecast projection time of nn hrs

Water temperature at forecast projection time of nn hrs

Note: The uppercase letters indicate a prediction series (e.g., H), and lowercase letters (e.g., h) indicate a reference series (observation or astro-nomical prediction). Slack water is defined as a current speed less than 0.26 m/s. The direction is computed only for current speeds greater than 0.26 m/s.

the forecast is valid. If there are, for example, 2 forecasts per day, then there will be two 6 h projec-tion values, separated by 12 hours in time. 2.3 Skill assessment statistics

Although no single set of statistics can quantify model performance perfectly, we have chosen sev-eral, easily-calculated quantities that provide rele-vant information on important categories of model behavior. A summary of relevant terms is shown in Table 2. For global assessment of errors, both the series mean (SM) and the frequency with which er-rors lie within specified limits (herein termed the central frequency, CF) are used. The SM indicates how well the model reproduces the observed mean and the CF indicates how often the error is within acceptable limits. The root mean square error (RMSE) and standard deviation (SD) are calculated, but have limited use since we do not expect errors to be normally distributed and CF is easier to explain to users lacking technical background. The CF con-cept has been used previously in National Ocean Service (NOS) of National Oceanic and Atmos-pheric Administration (NOAA) for data quality as-surance standards (Williams et al., 1989). The fre-quency of times of poor performance is determined by analyzing the outliers, which are values that ex-ceed specified limits. The positive outlier frequency (POF) measures how often the nowcast/forecast is higher than the observed. The maximum duration of positive outliers (MDPO) indicates whether there are long periods when the model overpredicts. The Negative Outlier Frequency (NOF) measures how often the nowcast/forecast is lower than the ob-served. The maximum duration of negative outliers (MDNO) indicates whether there are long periods when the model underpredicts. The MDPO and MDNO will be computed with data without gaps. For water levels, the “worst case”, from a model-based nowcast/forecast viewpoint, is when actual water level turns out to be low but the model erroneously predicted much higher water levels and the user would have been better off using the astro-nomical tide water level prediction. This is called the worst case outlier frequency (WOF).

30CHIN. J. OCEANOL. LIMNOL., 25(1), 2007 Vol.25

Table 2 Skill Assessment Statistics

Variable Explanation

The error is defined as the predicted value, p, minus the reference (observed or astronomical tide value, r: ei = pi

Error

- ri.

1N

Series Mean. The mean value of a series y. Calculated as y=∑yi SM Ni=1RMSE

Root Mean Square Error. Calculated as RMSE=

SD

CF(X) POF(X) NOF(X) MDPO(X)

Standard Deviation. Calculated as SD

Central Frequency. Fraction (percentage) of errors that lie within the limits ±X. Positive Outlier Frequency. Fraction (percentage) of errors that are greater than X. Negative Outlier Frequency. Fraction (percentage) of errors that are less than -X.

Maximum Duration of Positive Outliers. A positive outlier event is two or more consecutive occurrences of an error greater than X. MDPO is the length of time (based on the number of consecutive occurrences) of the long-est event.

Maximum Duration of Negative Outliers. A negative outlier event is two or more consecutive occurrences of an error less than -X. MDNO is the length of time (based on the number of consecutive occurrences) of the longest event.

Worst Case Outlier Frequency. Fraction (percentage) of errors that, given an error of magnitude exceeding X, either (1) the simulated value of water level is greater than the astronomical tide and the observed value is less than the astronomical tide, or (2) the simulated value of water level is less than the astronomical tide and the observed value is greater than the astronomical tide.

MDNO(X)

WOF(X)

2.4 Target frequency and duration limit

Most of the statistics have an associated target frequency of occurrence. For example,

S(X)≤P

where S is the statistics; X is the acceptable error magnitude (defined by the user); and P is the target frequency (or percentage).

CF(X)≥90%, POF(2X)≤1%, NOF(2X)≤1%

Other statistics are expressed as limits on the du-ration of errors, such as

MDPO(2X)≤L, MDNO(2X)≤L, WOF(2X)≤0.5%

where L is the time limit or maximum allowable duration.

According to the technical report (Hess, 2003) of the NOS in the NOAA, the values of acceptance error limits (X) and duration limits of the variables defined in Table 1 are summarized in Table 3, and these values are used for evaluating NOS’s opera-

tional nowcast and forecast hydrodynamic model systems.

3 DATA ANALYSIS TECHNIQUES

Observational and modeled data are processed and analyzed using several techniques. Observed time series that have gaps in the data are filled in one of three possible ways. Model-generated series, which are usually produced from numerous indi-vidual runs, must be concatenated to form a con-tinuous series. For each type of data, the entire se-ries is analyzed for harmonic constants and extrema (e.g., high water, maximum flood current) values. Specific methods for each process are discussed be-low.

Table 3 Acceptance error limit (X) and duration limit (L)

Variable* H, Hnn, AHW, ALW THW, TLW

U, Unn, AFC, AEC TFC, TEC

TSF, TEF, TSE, TEE D, Dnn DFC, DEC T, Tnn S, Snn

*

X L (hours) 15 cm 24 0.5 hours 25 0.26 m/s 24 0.5 hours 25 0.25 hours 25 22.5 degrees 24 22.5 degrees 25

7.7 degree ℃ 24 3.5 psu 24

Please refer to Table 1.

No.1 ZHANG et al.: The standards for skill assessment of operational marine forecast system 31

3.1 Gap filling and time interval conversion Data gaps often exist in observations. The ex-traction of extrema cannot be accomplished in a time series with gaps. Data gaps can be filled using different interpolation method. Three methods, (lin-ear interpolation, cubic spline interpolation, and singular value decomposition (SVD)) are adopted in the gap-filling program. As an option, the user can choose any method according to his/her experience and data simulation. If a gap is small enough, simple linear interpolation is appropriate. If a gap is large, a cubic spline or SVD interpolation should be used. The cubic spline interpolation is smooth in the first derivative, and continuous in the second derivative, both within an interval and at its boundaries. SVD produces a solution that is the best approximation in the least-squares sense in the case of an over-determined system (i.e., where the number of data points is greater than number of parameters), and SVD also produces a solution whose values are smallest in the least-squares sense in the case of an underdetermined system (i.e., where the number of data points is less than number of parameters, or if ambiguous combinations of parameters exist). SVD’s disadvantage is that it requires more memory space and can be significantly slower than solving the normal equations. However, its great advantage is that it (theoretically) cannot fail, and this more than makes up for the speed disadvantage.

The time intervals of observation and modeled time series might be different, so all time series with different time intervals have to be converted into equally-spaced time series with the same unique de-sired time interval. 3.2 Filtering

Because of short period variations and noise, filtering of values in a time series is sometimes nec-essary to select accurately the extrema (i.e., maxi-mum and minimum) values and times. A Fourier filter is a good choice as it computes the amplitudes of the components of the signal at various frequen-cies and reduces the amplitudes at selected frequen-cies. Simple smoothing should be avoided because it reduces extrema amplitudes.

3.3 Tidal prediction and harmonic analysis

Tidal prediction of water level and current is

required for skill assessment in tidal regions. Tidal harmonic constants can be obtained either from some government agencies or can be derived from observations or model outputs using Fourier har-monic analysis or least squares harmonic analysis methods. Astronomical tidal water level and current time series will be predicted from 37 tidal constitu-ents for any time period. The 37 constituents used by in this paper are M2, S2, N2, K1, M4, O1, M6, MK3, S4, MN4, ν2, S6, µ2, 2N2, OO1, λ2, S1, M1, J1, Mm, Ssa, Sa, Msf, Mf, ρ1, Q1, T2, R2, 2Q1, P1, 2SM2, M3, L2, 2MK3, K2, M8, and MS4. A description of the constituents can be found in Schureman (1958).

In tidal regions, a comparison of tidal harmonic constants is necessary for the evaluation of water levels and currents. For this comparison, harmonic constants (amplitudes and phases) are analyzed from tide-only model simulation and observed data. Two analytical techniques, least squares harmonic analy-sis and Fourier harmonic analysis, are used in terms of the length of the data time series. The least squares method (Zervas, 1999) is a method for de-riving the tidal constituents from a water level or current time series by creating a matrix of covari-ance between each individual constituent time series and the observed time series. The matrix is inverted to solve for the amplitudes and phases of the har-monic constituents. The constituent with the highest correlation is then subtracted from the observed time series, and the matrix is recalculated with a re-sidual time series in place of the observed. This method can solve 175 tidal constituents, but will not analyze less than 29 days of data. The Fourier har-monic analysis method (Dennis and Long, 1971) uses Fourier series summations to obtain the tidal constituents of water level or current data. This method has been programed for data periods of ei-ther 15 or 29 days of continuous data time series. 3.4 Extrema extraction

For skill assessment, the amplitudes and times of high/low waters, and those of maximum flood/ebb currents are required. The time series needs to be filtered if there is noise before extracting extrema. The extrema are extracted by searching for the largest and smallest values within a given time period in a series by the following method. First, the time series values within each 30 min segment are averaged to obtain a new series with a time interval