XU Jiangling, HUANG Juan, GAO Song, and CAO Yajing

North China Sea Marine Forecast Center of the State Oceanic Administration, Qingdao 266061, P. R. China

Assimilation of High Frequency Radar Data into a Shelf Sea Circulation Model

XU Jiangling*, HUANG Juan, GAO Song, and CAO Yajing

North China Sea Marine Forecast Center of the State Oceanic Administration, Qingdao 266061, P. R. China

High Frequency (HF) radar current data is assimilated into a shelf sea circulation model based on optimal interpolation (OI) method. The purpose of this work is to develop a real-time computationally highly efficient assimilation method to improve the forecast of shelf current. Since the true state of the ocean is not known, the specification of background error covariance is arduous. Usually, it is assumed or calculated from an ensemble of model states and is kept in constant. In our method, the spatial covariances of model forecast errors are derived from differences between the adjacent model forecast fields, which serve as the forecast tendencies. The assumption behind this is that forecast errors can resemble forecast tendencies, since variances are large when fields change quickly and small when fields change slowly. The implementation of HF radar data assimilation is found to yield good information for analyses. After assimilation, the root-mean-square error of model decreases significantly. Besides, three assimilation runs with variational observation density are implemented. The comparison of them indicates that the pattern described by observations is much more important than the amount of observations. It is more useful to expand the scope of observations than to increase the spatial interval. From our tests, the spatial interval of observation can be 5 times bigger than that of model grid.

data assimilation; current radar; shelf circulation model

1 Introduction

With the rapid development of HF radar (Barrick, 1978; Barricket al., 1977), the coastal surface current can be observed more reliably with a high temporal and spatial resolution, compared toin-situdata, even to satellite altimetry data. For this reason more and more assimilation efforts have been made in coastal current studies.

The earliest work in this context was performed by Lewiset al.(1998), who treated the Doppler radar surface current data as an additional layer of moving water overlying the ocean surface, and used the differences between the observed values and the model velocity to correct the wind forcing data for the model. By this approach the main pattern of the Doppler radar currents was introduced into the model surface current calculation.

Scottet al.(2000) used an idealized, linear model with a variational inverse assimilation scheme to assess the feasibility of assimilation of surface current measurements. Kurapovet al.(2002) did a comparison of the generalized inverse method (GIM), Kalman Filter (KF) and optimal interpolation (OI) methods with an idealized three-dimensional time-dependent coastal baroclinic model. They found that the GIM yielded a smaller posterior error variance than the KF or OI method because the information of the previous flow fields is introduced into the data. Compared to Lewiset al.(1998), the assimilation methods described in these two articles are more sophisticated; however, the models are still linear,i.e., not fully realistic.

In order to perform realistic simulations, Breivik and Saetra (2001) assimilated HF radar-derived current data into a realistic coastal model using an adjusted optimal interpolation method, in which an ensemble of model states was sampled from a reference run in order to calculate the required background error covariance matrix. They showed that a significant improvement on the model forecast by assimilation of HF radar currents can last for six hours.

Okeet al.(2002) used a sequential optimal interpolation scheme to assimilate HF radar data into a primitive equation coastal ocean model. An ensemble of model states from 18 different summers was used to compute the forecast error covariances, and a time-averaging procedure was used to introduce the corrections gradually to the model state. The same assimilation method was applied by Paduan and Shulman (2004) to assimilate HF radar data observed in the Monterey Bay area. The assimilation analysis was only implemented for surface variables; by an additional consideration of an Ekmanlayer projection the spatial and temporal differences between modeled and observed data decreased significantly.

More recently, Barthet al.(2008) developed an en-semble method. They used a nested model under different wind forcing conditions to estimate the error covariance of the model state vectors and the covariance between ocean currents and winds. In this study tides were removed from the observed surface currents, which, moreover, were averaged over two days. Chenet al.(2011) implemented a Quick-EnKF method in the Changjiang Estuary. That assimilation method is essentially equal to OI method, while with a variable background error covariance.

All of these studies demonstrate that OI method is a nice choice for practical reason. Otherwise, computational costs can become extremely high. Besides, OI method has been used in many different ways, like those in Krishnamurthyet al.(2010), Loyola and Coldewey-Egbers (2012) and so on. However, the way to get an estimation of model error covariance used in OI is varying,e.g., an ensemble of model states for different summers or only different wind forcing conditions are employed. In this paper, the difference between the adjacent background fields serving as the forecast tendency is used to get perturbation matrix, and then the analysis field is calculated in such cases that the covariance between observational error and perturbation errors comes to a minimum. The general assumption behind this implementation is the same as that for the Canadian Quick Covariance (CQC, (Polavarapuet al., 2005; Jacksonet al., 2008)) method, namely, the forecast tendencies can resemble forecast errors. It has to be noted that the main objective of the current study was to develop a computationally efficient system to assimilate HF radar data into a quasioperational shelf sea circulation model. The outline of this paper is as follows: The assimilation system is briefly described in Section 2, with the information concerning the radar data and the applied model also given herein. In Section 3 the design of assimilation experiments and results are presented. The summary is given in Section 4.

2 Data Assimilation System

2.1 Observation: HF Radar Data

The HF radar system called CODAR is mounted on the islands of Xuejiadao and Xiaomaidao off the east coast of Qingdao. Each radar array measures the radial component of the surface current. Each radial current velocity is a statistical result from samples of velocities derived from lines in the Doppler spectrum (Gurgel and Antonischki, 1997; Gurgelet al., 1999). CODAR measures surface currents within an area of about 20 km by 20 km with a spatial resolution of 0.5 km and with temporal sampling of one hour. The radar systems are deployed along the coast and measure radial current speeds. Hence, two sites are necessary for composing a two-dimensional current vector. The surface current measured is a horizontal mean over several km2(range and azimuthal resolution), over about the upper 0.5–1 m of the water.

The data we used are already converted in eastward and northward directions. Fig.1 gives an overview on the area and typical radar coverage during the experiment. The radar coverage changed because of radio interference, sun spot activity, varying sea state,etc. Besides, the geometric error is generated when combing two radial current components, and sequentially has been translated tou- andv- components of the radar current vectors. Since the accuracy is calculated from the variance of the radial current velocity, it should be the diagonal part of error covariance matrix of the measurements.

Fig.1 An overview of the typical radar coverage measured by the CODAR.

2.2 Model Configuration

HAMburg Shelf Ocean Model (HAMSOM) is a high resolution 3-D baroclinic, free surface, shallow water equation model. Its numerical scheme is described by Backhaus (1985), with further modification in the turbulent closure scheme by Pohlmann (1996a, b). The horizontal diffusion is determined by the Smagorinsky scheme (Smagorinsky, 1963). It is known that assimilation analysis consumes much CPU time, so one of the reason we chose HAMSOM is because of the parallelization of this code, which in this model is performed bymeans of the domain splitting method (Pohlmann, 2006). Meanwhile we also parallelized the assimilation procedure.

In the present study, the model horizontal resolution is 1/240 degree (see Fig.2). Thez-coordinate system is applied in the vertical direction using a maximum of 6 layers. The resolution in the upper 12 m is 4 m, and that in the last three layers is 8 m. The model time step is 30 second. The model surface forcing parameters, such as heat flux, air temperature, relative humidity, cloud cover, precipitation, as well as wind are supplied by NCEP (Kalnayet al., 1996). All these variables are given every 6 hours, and the initial temperature and salinity fields are obtained from the climatological LEVITUS data set (Levituset al., 1994). As known, tides are dominating elements for the shelf circulation in our model domain, so the tidal influence is definitely included in our model simulation by adding the tidal elevation to the sea level height at the open boundaries. The tidal elevation is calculated according to amplitudes and phases of M2, S2, K1, O1, which are extracted from OTPS (OSU Tidal Prediction Software).

Fig.2 Model domain topography (m).

2.3 Optimal Interpolation

Optimal Interpolation is a sequential data assimilation method using predefined error statistics, which are time invariant. The OI is an algebraic simplification of the computation of the weightKin the analysis equations given as

here,Xis the n-dimensional model state vector, with subscripts ‘a(chǎn)’ and ‘b’ denoting analysis and background;Kis the gain matrix, which is ratio between background (numerical model) error covariance (B) and observation error covariance (R), andyis the observational data vector;His the observation operator, which maps the model state onto the observation space; T donates the transposition of matrix.

In the OI algorithm it is necessary to have the background error covariances so that can easily be applied to pairs of model and observed variables, or to pairs of observed variables. However, it is difficult to implement if the observation operator is complex. Moreover, since the true state will never be known, the specification of background error covariance is arduous. In fact, even the more sophisticated methods, like EnKF and variation methods, are established on the treatment ofB(Buehneret al., 2010; Law and Stuart, 2012). All the previous efforts related to radar data assimilation are already described in Section 1. In our quick assimilation method the model is integrated only once and thus the difference between the adjacent background fields serves as the forecast tendency to get perturbation matrix. For example, if the number of the ensemble is set toN, in the process of integrationNdifferences between two adjacent backgroundare used to get the background error covariance matrixThis variable helps to calculate the weighting of background data. The assumption of this usage agrees with the Canadian Quick Covariance method which regards the tendency of an integrated process as a forecast error. More details were described by (Chen and Zhan, 2011). Then the analysis field is calculated in such cases that the covariance between observational error and perturbation error comes to a minimum.

3 Assimilation Implementation and Results

Twins experiments are performed: one of them is without assimilation and called ‘Freerun’, the other is performed with our OI method and named as Assim1. As mentioned before, the spatial resolution of the radar data is 500 m, which is comparable to our model spatial resolution.

All these model simulations are initialized on January 1, 2012, and the assimilation procedure starts from February 1 till end of the month. The analysis time interval is generally one hour. The reason we used the word ‘generally’is because sometimes the number of observations is very small, though that situation occurs not very often. Only the variables at the surface are analyzed, since the radar data only represents the very thin surface layer, and the vertical correlations between these variables are poor (pictures are not shown here). Besides, all the variables of HAMSOM are calculated from the surface to the bottom. Hence, after the surface variables are updated, the vertical projection of the update can be done by the propagation of model itself.

3.1 Analysis of Assimilated Performance

First and foremost, the comparison of surface current fields without (Freerun) and with assimilation (Assim1) at one random moment is presented in Fig.3 to evaluate the influence of radar data assimilation. The figure shows that the primary updates by assimilation exist in the radar domain and its vicinity. However, with the integration ofnumerical model, currents in the whole model domain are modified. Compared to the radar data (right panel in Fig.1), we found that our free model results do not agree well with the radar data, since the main direction of velocity of radar data is eastwards, while that of Freerun in the same area is northeastwards. Besides, the magnitude of velocity of radar data is a little bigger than that of the reference run. From this picture, we can also see that our new assimilation method gives an obvious improvement in both velocity direction and magnitude. It is obvious that the assimilation scheme captures the strength and extent of the coastal current very well at this particular instant.

In order to check the added information by assimilation, we computed the spatially averaged kinetic energy in the radar-covered area for both model and observations:

where ? is the area covered by the radar andUis the horizontal current speed. The density is assumed constant and left out. This method has the advantage of appreciating a corrected current field (typically the coastal current) even when the current maximum is slightly dislocated by the model yet still improved over the free run.

Fig.3 Comparison of surface current fields without (blue) and with (red) assimilation at one random moment (02:00 on February 15).

Fig.4 Scatter plots of observed and modeled current speeds from the free run (upper panel) and the run with assimilation (lower panel). The correlation coefficient r is given in the upper left corner of each panel. The linear regression best line is indicated together with the ideal 45? line.

The scatter plot is a more conventional way to compare datasets. Fig.4 shows scatter plots of observed and modeled current speeds from the assimilation (lower panel) and the free run (upper panel). The correlation coefficient‘r’ is given in the upper left-hand corner of each panel. The linear regression best-fit line is indicated together with the ideal 45? line. As can be seen, the free run underestimates the energy level of the coastal current. Assimilation results in a significant improvement over the free run, with an energy level on a par with what is observed.In order to assess the quality of assimilation with our OI method over time, the root-mean-square (RMS) errors of bothUandVcomponents in Assim1 and the reference run relative to the radar data are computed in the observational space (Fig.5). The pictures indicate that the RMS errors of both components in the assimilation run are much smaller than in the reference run. After assimilation, the errors are reduced to around 10 cm s-1. Besides, the comparison of RMS errors of theU-velocity andV- velocity indicates that the correction of the former variable is smaller than that of the latter variable. This can be explained by the fact that in the radar domain, the alongcoast velocity is bigger than the cross-coast one.

In a word, assimilation using radar data makes the surface velocity structure agree better with the radar observations.

Fig.5 Root-mean-square errors ofU-component (upper panel) andV-component (lower panel) relative to radar data, in the reference run (blue line) and Assim1 (red line), respectively.

3.2 Variation of Updated Surface Currents

As mentioned before, the spatial resolution of the radar data is comparable to our model spatial resolution. However, normally the spatial interval of observation from buoy or satellite is big, and we wonder what spatial resolution of observation can be appropriate for data assimilation. Hence, with this variation of the observational density, another two experiments were conducted, and compared to the reference run. The details of these experiments are described as follows:

Freerun: The reference run, which is already described before.

Assim1: Assimilation run. In this experiment, all of the available observational data are used. The cutoff radius is set to 0.03?, which is almost six times the size of model grid interval.

Assim2: The same assimilation method and setup used as in Assim1 run but with only 20% observations.

Assim3: The same assimilation method and setup as in Assim1 run but with only 5% observations.

First, the comparison of velocity fields with the variation of the observational density is presented. Fig.6 shows surface current fields zoomed into the radar area of these two experiments at the same moment as in Assm1. When comparing this figure with Fig.3, only a slight difference can be found. It seems that the velocity corrections by the assimilation decrease with the number of observations decreasing. Even with 20% observational data, the simulated velocity field can still be improved. However, when only 50% of observations are used, the effect of assimilation decreases significantly, and the background field (result of reference run) plays an increasing role in the assimilation. To make a quantitative analysis, RMS errors of these four experiments at this moment are calculated and listed in Table 1.

Table 1 RMS errors (cm s-1)

To assess the changes caused by the assimilation over time, we calculated the RMS errors and energy of these three assimilation experiments relative to the HF radar data, which are shown in Fig.7 and Fig.8 respectively. The pictures clearly show that as the amount of observations decreases, the RMS errors increase, and the assimilation effects decrease. It is also obvious that there are big differences between assimilation run with all observations and assimilation run with only 5% observations, but there are small differences between Assim1 and Assim2. It means that the pattern described by observations is more important than the amount of observations.

Fig.6 Comparison of surface current fields without (blue) and with (red) assimilation at one random moment (02:00 on February 15) of Assim2 (upper panel) and Assim3 (lower panel), respectively.

Fig.7 Root-mean-square errors of U-component (upper panel) and V-component (lower panel) relative to radar data in the reference run (blue), Assim1 (red) , Assim2 (pink) and Assim3 (green), respectively.

Fig.8 Scatter plots of observed and modeled current speeds from the Assim2 (upper panel) and Assim3 run (lower panel). The correlation coefficient r is given in the upper left corner of each panel. The linear regression best line is indicated together with the ideal 45? line.

4 Conclusions

In the present study, HF radar current data is assimilated into a shelf sea circulation model based on optimal interpolation method. The implementation of modeling shows that this method is suitable for operational assimilation, and makes great improvement in shelf current simulation.

The comparisons between the three assimilation runs with different observational densities indicate that, firstly, our OI assimilation method is successful, since with the number of observations decreasing, the effect of assimilation decreases too. That can be easily deduced from the changes of spatial RMS errors. Secondly, the pattern described by observations is more important than the amount of observations. It would be more useful to expand the scope of observations than to increase the spatial interval.

Acknowledgements

This work is supported by the State Oceanic Administration Young Marine Science Foundation (No. 2013201), the Shandong Provincial Key Laboratory of Marine Ecology and Environment & Disaster Prevention and Mitigation Foundation (No. 2012007), the Marine Public Foundation (No. 201005018) and the North China Sea Branch Scientific Foundation (No. 2014B10). The authors are grateful to Dr. Lingjuan Wu for her suggestions. Many thanks also go to the reviewers for their helpful comments.

Backhaus, J., 1985. A three-dimensional model for the simulation of shelf sea dynamics. Ocean Dynamics, 38 (4): 165-187.

Barrick, D., 1978. HF radio oceanography review. Boundary-Layer Meteorology, 13 (1): 23-43.

Barrick, D., Evans, M., and Weber, B., 1977. Ocean surface currents mapped by radar. Science, 198 (4313): 138-144.

Barth, A., Alvera-Azcarate, A., and Weisberg, R., 2008. Assimilation of high-frequency radar currents in a nested model of the west Florida shelf. Journal of Geophysical Research, 113 (C8), DOI: 10.1029/2007JC004585.

Breivik, R., and Saetra, R., 2001. Real time assimilation of HF radar currents into a coastal ocean model. Journal of Marine Systems, 28 (3-4): 161-182.

Buehner, M., Houtekamer, P. L., Charette, C., and Mitchell, H. L., 2010. Intercomparison of variational data assimilation and the ensemble Kalman filter for global deterministic NWP. Part I: Description and single-observation experiments. Monthly Weather Review, 138 (5): 1550-1566.

Chen, X. E., and Zhan, P., Chen, J. R., and Qian, H. B., 2011. Numerical study of current fields near the Changjiang Estuary and impact of Quick-EnKF assimilation. Acta Oceanologica Sinica, 30 (5): 33-44.

Gurgel, K., and Antonischki, G., 1997. Measurement of surface current fields with high spatial resolution by the HF radar WERA. Geoscience and Remote Sensing, 1997. IGARSS’97.

Remote Sensing–A Scientific Vision for Sustainable Development, 1997 IEEE International Vol. 4, 1820-1822.

Gurgel, K., Antonischki, G., Essen, H., and Schlick, T., 1999. Wellen Radar (WERA): A new ground-wave HF radar for ocean remote sensing. Coastal Engineering, 37 (3-4): 219-234.

Jackson, D. R., Keil, M., and Devenish, B. J., 2008. Use of Canadian quick covariances in the met office data assimilation system. Quarterly Journal of Royal Meteorological Society, 134: 1567-1582.

Kalnay, E., Kanamitsu, M., Kistler, R., Collins, W., Deaven, D., Gandin, L., Iredell, M., Saha, S., White, G., Woollen, J., and Zhu, Y., 1996. The NCEP/NCAR 40-year reanalysis project. Bulletin of the American Meteorological Society, 77 (3): 437-471.

Krishnamurthy, A., Cobb, L., Mandel, J., and Beezley, J., 2010. Bayesian tracking of emerging epidemics using ensemble optimal statistical interpolation. 2010 Joint Statistical Meetings, July 31–August 5, 2010, Vancouver, Canada.

Kurapov, A., Egbert, G., Miller, R., and Allen, J., 2002. Data assimilation in a baroclinic coastal ocean model: Ensemble statistics and comparison of methods. Monthly Weather Review, 130 (4): 1009-1025.

Law, K. J. H., and Stuart, A. M., 2012. Evaluating data assimilation algorithms. Monthly Weather Review, 140: 3757-3782.

Levitus, S., Boyer, T., and Antonov, J., 1994. World Ocean Atlas 1994. Volume 5. Interannual Variability of Upper Ocean Thermal Structure. Technique report, PB-95-270120/ XAB, National Environmental Satellite, Data, and Information Service, Washington, DC, 176pp.

Lewis, J., Shulman, I., and Blumberg, A., 1998. Assimilation of CODAR observations into ocean models. Continental Shelf Research, 18: 541-559.

Loyola, D. G., and Coldewey-Egbers, M., 2012. Multi-sensor data merging with stacked neural networks for creation of satellite long-term climate data records. EURASIP Journal on Advances in Signal Processing, 2012: 91.

Oke, P., Allen, J., Miller, R., Egbert, G., and Kosro, P., 2002. Assimilation of surface velocity data into a primitive equation coastal ocean model. Journal of Geophysical Research, 107, 3122, DOI: 10.1029/2000JC000511.

Paduan, J., and Shulman, I., 2004. HF radar data assimilation in the Monterey Bay area. Journal of Geophysical Research, 109 (C7), DOI: 10.1029/2003JC001949.

Pohlmann, T., 1996a. Calculating the development of the thermal vertical stratification in the North Sea with a three-dimensional baroclinic circulation model. Continental Shelf Research, 16 (2): 163-194.

Pohlmann, T., 1996b. Predicting the thermocline in a circulation model of the North seapart I: Model description, calibration and verification. Continental Shelf Research, 16 (2): 131-146.

Pohlmann, T., 2006. A meso-scale model of the central and southern North Sea: Consequences of an improved resolution. Continental Shelf Research, 26 (19): 2367-2385.

Polavarapu, S., Ren, S., Rochon, Y., Sankey, D., Ek, N., Koshyk, J., and Tarasick, D., 2005. Data assimilation with the Canadian middle atmosphere model. Atmosphere-Ocean, 43 (1): 77-100.

Scott, R., Allen, J., Egbert, G., and Miller, R., 2000. Assimilation of surface current measurements in a coastal ocean model. Journal of Physical Oceanography, 30 (9): 2359-2378.

Smagorinsky, J., 1963. General circulation experiments with the primitive equations. I. The basic experiment. Monthly Weather Review, 91: 99-164.

(Edited by Xie Jun)

(Received December 13, 2012; revised March 18, 2013; accepted December 6, 2013)

? Ocean University of China, Science Press and Springer-Verlag Berlin Heidelberg 2014

* Corresponding author. Tel: 0086-532-58750655

E-mail: xujiangling@nmfc.gov.cn