Data assimilation with a barotropically unstable shallow water system using representer algorithms

The cycling representer algorithm of Xu and Daley (2000) is a weak constraint four-dimensional variational data assimilation algorithm. It was successfully applied to a one-dimensional transport problem and was able to successfully extract the signal from noisy and sparse observations. The algorithm, however, has not previously been applied to a multivariate, multidimensional system with dynamic instability. The algorithm is also very computationally demanding and awaits considerable enhancement in computer power before being practical for operational forecast models. We have two objectives in this paper. The first is to apply the cycling representer algorithm to a two-dimensional, multivariate barotropically unstable linear shallow water system. The second objective is to formulate and test an accelerated representer algorithm that is much more computationally tractable than the cycling representer algorithm itself. A linear shallow water system with a barotropically unstable basic state was used as a test bed to conduct data assimilation experiments. The evolution of a 'neutral' eastward-propagating singular vector was selected as the 'truth', against which all data assimilation experiments were to be evaluated. The results indicated that the cycling representer algorithm was capable of providing satisfying state estimates for a multivariate, multidimensional system. The results from the accelerated representer algorithm were very encouraging because it is sufficiently computationally tractable to be used on present day multi-processor machines for operational applications.     

Time‐space weak‐constraint data assimilation for nonlinear models

A general perturbation–linearization scheme is proposed for the problem of data assimilation with an imperfect and nonlinear model, allowing for the application of the weak constraint representer method. The scheme is shown in discrete formalism for a generic model. An application example is given with computer‐generated data in the case of the Burgers equation. Discussion in reference to the assimilation example concerns: the rôle of the model error, seen as a forcing term in the dynamics; the rôle of representers as a posteriori error covariances; a comparison among different choices for a priori dynamic error variance and strong constraint assimilation. Weak and strong constraint methods are also compared in a forecasting experiment.     

Cycling the representer algorithm for variational data assimilation with a nonlinear reduced gravity ocean model

The representer method was used by [Ngodock, H.E., Jacobs, G.A., Chen, M., 2006. The representer method, the ensemble Kalman filter and the ensemble Kalman smoother: a comparison study using a nonlinear reduced gravity ocean model. Ocean Modelling 12, 378–400] in a comparison study with the ensemble Kalman filter and smoother involving a 1.5 nonlinear reduced gravity idealized ocean model simulating the Loop Current (LC) and the Loop Current eddies (LCE) in the Gulf of Mexico. It was reported that the representer method was more accurate than its ensemble counterparts, yet it had difficulties fitting the data in the last month of the 4-month assimilation window when the data density was significantly decreased. The authors attributed this failure to increased advective nonlinearities in the presence of an eddy shedding causing the tangent linear model (TLM) to become inaccurate. In a separate study [Ngodock, H.E., Smith, S.R., Jacobs, G.A., 2007. Cycling the representer algorithm for variational data assimilation with the Lorenz attractor. Monthly Weather Review 135 (2), 373–386] applied the cycling representer algorithm to the Lorenz attractor and demonstrated that the cycling solution was able to accurately fit the data within each cycle and beyond the range of accuracy of the TLM, once adjustments were made in the early cycles, thus overcoming the difficulties of the non-cycling solution. The cycling algorithm is used here in assimilation experiments with the nonlinear reduced gravity model. It is shown that the cycling solution overcomes the difficulties encountered by the non-cycling solution due to a limited time range of accuracy of the TLM. Thus, for variational assimilation applications where the TLM accuracy is limited in time, the cycling representer becomes a very powerful and attractive alternative, given that its computational cost is significantly lower than that of the non-cycling algorithm.     

Cycling the representer algorithm for variational data assimilation with a nonlinear reduced gravity ocean model

The representer method was used by [Ngodock, H.E., Jacobs, G.A., Chen, M., 2006. The representer method, the ensemble Kalman filter and the ensemble Kalman smoother: a comparison study using a nonlinear reduced gravity ocean model. Ocean Modelling 12, 378–400] in a comparison study with the ensemble Kalman filter and smoother involving a 1.5 nonlinear reduced gravity idealized ocean model simulating the Loop Current (LC) and the Loop Current eddies (LCE) in the Gulf of Mexico. It was reported that the representer method was more accurate than its ensemble counterparts, yet it had difficulties fitting the data in the last month of the 4-month assimilation window when the data density was significantly decreased. The authors attributed this failure to increased advective nonlinearities in the presence of an eddy shedding causing the tangent linear model (TLM) to become inaccurate. In a separate study [Ngodock, H.E., Smith, S.R., Jacobs, G.A., 2007. Cycling the representer algorithm for variational data assimilation with the Lorenz attractor. Monthly Weather Review 135 (2), 373–386] applied the cycling representer algorithm to the Lorenz attractor and demonstrated that the cycling solution was able to accurately fit the data within each cycle and beyond the range of accuracy of the TLM, once adjustments were made in the early cycles, thus overcoming the difficulties of the non-cycling solution. The cycling algorithm is used here in assimilation experiments with the nonlinear reduced gravity model. It is shown that the cycling solution overcomes the difficulties encountered by the non-cycling solution due to a limited time range of accuracy of the TLM. Thus, for variational assimilation applications where the TLM accuracy is limited in time, the cycling representer becomes a very powerful and attractive alternative, given that its computational cost is significantly lower than that of the non-cycling algorithm.     

The representer method, the ensemble Kalman filter and the ensemble Kalman smoother: A comparison study using a nonlinear reduced gravity ocean model

This paper compares contending advanced data assimilation algorithms using the same dynamical model and measurements. Assimilation experiments use the ensemble Kalman filter (EnKF), the ensemble Kalman smoother (EnKS) and the representer method involving a nonlinear model and synthetic measurements of a mesoscale eddy. Twin model experiments provide the “truth” and assimilated state. The difference between truth and assimilation state is a mispositioning of an eddy in the initial state affected by a temporal shift. The systems are constructed to represent the dynamics, error covariances and data density as similarly as possible, though because of the differing assumptions in the system derivations subtle differences do occur. The results reflect some of these differences in the tangent linear assumption made in the representer adjoint and the temporal covariance of the EnKF, which does not correct initial condition errors. These differences are assessed through the accuracy of each method as a function of measurement density. Results indicate that these methods are comparably accurate for sufficiently dense measurement networks; and each is able to correct the position of a purposefully misplaced mesoscale eddy. As measurement density is decreased, the EnKS and the representer method retain accuracy longer than the EnKF. While the representer method is more accurate than the sequential methods within the time period covered by the observations (particularly during the first part of the assimilation time), the representer method is less accurate during later times and during the forecast time period for sparse networks as the tangent linear assumption becomes less accurate. Furthermore, the representer method proves to be significantly more costly (2–4 times) than the EnKS and EnKF even with only a few outer iterations of the iterated indirect representer method.     

Weak and strong constraint data assimilation in the inverse Regional Ocean Modeling System (ROMS): Development and application for a baroclinic coastal upwelling system

We describe the development and preliminary application of the inverse Regional Ocean Modeling System (ROMS), a four dimensional variational (4DVAR) data assimilation system for high-resolution basin-wide and coastal oceanic flows. Inverse ROMS makes use of the recently developed perturbation tangent linear (TL), representer tangent linear (RP) and adjoint (AD) models to implement an indirect representer-based generalized inverse modeling system. This modeling framework is modular. The TL, RP and AD models are used as stand-alone sub-models within the Inverse Ocean Modeling (IOM) system described in [Chua, B.S., Bennett, A.F., 2001. An inverse ocean modeling system. Ocean Modell. 3, 137–165.]. The system allows the assimilation of a wide range of observation types and uses an iterative algorithm to solve nonlinear assimilation problems. The assimilation is performed either under the perfect model assumption (strong constraint) or by also allowing for errors in the model dynamics (weak constraints). For the weak constraint case the TL and RP models are modified to include additional forcing terms on the right hand side of the model equations. These terms are needed to account for errors in the model dynamics.Inverse ROMS is tested in a realistic 3D baroclinic upwelling system with complex bottom topography, characterized by strong mesoscale eddy variability. We assimilate synthetic data for upper ocean (0–450 m) temperatures and currents over a period of 10 days using both a high resolution and a spatially and temporally aliased sampling array. During the assimilation period the flow field undergoes substantial changes from the initial state. This allows the inverse solution to extract the dynamically active information from the synthetic observations and improve the trajectory of the model state beyond the assimilation window. Both the strong and weak constraint assimilation experiments show forecast skill greater than persistence and climatology during the 10–20 days after the last observation is assimilated.Further investigation in the functional form of the model error covariance and in the use of the representer tangent linear model may lead to improvement in the forecast skill.     

Reducing systematic errors by empirically correcting model errors

A methodology for the correction of systematic errors in a simplified atmospheric general‐circulation model is proposed. First, a method for estimating initial tendency model errors is developed, based on a 4‐dimensional variational assimilation of a long‐analysed dataset of observations in a simple quasi‐geostrophic baroclinic model. Then, a time variable potential vorticity source term is added as a forcing to the same model, in order to parameterize subgrid‐scale processes and unrepresented physical phenomena. This forcing term consists in a (large‐scale) flow dependent parametrization of the initial tendency model error computed by the variational assimilation. The flow dependency is given by an analogues technique which relies on the analysis dataset. Such empirical driving causes a substantial improvement of the model climatology, reducing its systematic error and improving its high frequency variability. Low‐frequency variability is also more realistic and the model shows a better reproduction of Euro‐Atlantic weather regimes. A link between the large‐scale flow and the model error is found only in the Euro‐Atlantic sector, other mechanisms being probably the origin of model error in other areas of the globe.     

Temporal evolution of innovation and residual statistics in the ECMWF variational data assimilation systems

The temporal evolution of innovation and residual statistics of the ECMWF 3D‐ and 4D‐Var data assimilation systems have been studied. First, the observational method is applied on an hourly basis to the innovation sequences in order to partition the perceived forecast error covariance into contributions from observation and background errors. The 4D‐Var background turns out to be ignificantly more accurate than the background in the 3D‐Var. The estimated forecast error variance associated with the 4D‐Var background trajectory increases over the assimilation window. There is also a marked broadening of the horizontal error covariance length scale over the assimilation window. Second, the standard deviation of the residuals, i.e., the fit of observations to the analysis is studied on an hourly basis over the assimilation window. This fit should, in theory, reveal the effect of model error in a strong constraint variational problem. A weakly convex curve is found for this fit implying that the perfect model assumption of 4D‐Var may be violated with as short an assimilation window as six hours. For improving the optimality of variational data assimilation systems, a sequence of retunes are needed, until the specified and diagnosed error covariances agree.     

Moisture retrievals from simulated zenith delay "observations" and their impact on short‐range precipitation forecasts

The feasibility of assimilating the GPS total zenith delay into atmospheric models is investigated within the framework of the "Observing System Simulation Experiment." The total zenith delay is made up of two terms: one is proportional to the pressure at the site of the GPS ground‐based receiver and the other to the overlying amount of water vapor. Using the MM5 mesoscale model and its adjoint, a set of 4‐dimensional variational (4DVAR) experiments is performed. Results from the assimilation of simulated precipitable water observations are used as the benchmark. The model domain covers Southern California. The observations are simulated with a 10 km horizontal resolution model that includes full physics, while a 20‐km resolution and a less comprehensive physics package are used in the 4DVAR experiments. Both, the 10‐km and 20‐km models employ the same set of 15 vertical levels. Moisture fields retrieved from the total zenith delay are found to compare very well with those retrieved from the precipitable water. Verified against the observations, the vertically integrated moisture is found to be very accurate. An overall improvement is also achieved in the vertical profiles of the moisture fields. The use of the so‐called background term and model initialization are shown to greatly reduce the negative impact that the sole assimilation of the total zenith delay can have on the pressure field and integrated water vapor. The adverse effect stems from the poor resolution of the topography needed to evaluate the model pressure at the GPS sites. The analysis increments of all model fields are found to be similar to the counterparts obtained from the assimilation of the precipitable water. The same is true for the short‐range precipitation forecasts initiated from the 4DVAR‐optimal initial conditions.     

Low‐dimensional representation of error covariance

Ensemble and reduced‐rank approaches to prediction and assimilation rely on low‐dimensional approximations of the estimation error covariances. Here stability properties of the forecast/analysis cycle for linear, time‐independent systems are used to identify factors that cause the steady‐state analysis error covariance to admit a low‐dimensional representation. A useful measure of forecast/analysis cycle stability is the bound matrix , a function of the dynamics, observation operator and assimilation method. Upper and lower estimates for the steady‐state analysis error covariance matrix eigenvalues are derived from the bound matrix. The estimates generalize to time‐dependent systems. If much of the steady‐state analysis error variance is due to a few dominant modes, the leading eigenvectors of the bound matrix approximate those of the steady‐state analysis error covariance matrix. The analytical results are illustrated in two numerical examples where the Kalman filter is carried to steady state. The first example uses the dynamics of a generalized advection equation exhibiting non‐modal transient growth. Failure to observe growing modes leads to increased steady‐state analysis error variances. Leading eigenvectors of the steady‐state analysis error covariance matrix are well approximated by leading eigenvectors of the bound matrix. The second example uses the dynamics of a damped baroclinic wave model. The leading eigenvectors of a lowest‐order approximation of the bound matrix are shown to approximate well the leading eigenvectors of the steady‐state analysis error covariance matrix.     

The Regional Ocean Modeling System (ROMS) 4-dimensional variational data assimilation systems: Part II - Performance and application to the California Current System

The Regional Ocean Modeling System (ROMS) 4-dimensional variational (4D-Var) data assimilation systems have been systematically applied to the mesoscale circulation environment of the California Current to demonstrate the performance and practical utility of the various components of ROMS 4D-Var. In particular, we present a comparison of three approaches to 4D-Var, namely: the primal formulation of the incremental strong constraint approach; the dual formulation “physical-space statistical analysis system”; and the dual formulation indirect representer approach. In agreement with theoretical considerations all three approaches converge to the same ocean circulation estimate when using the same observations and prior information. However, the rate of convergence of the dual formulation was found to be inferior to that of the primal formulation. Other aspects of the 4D-Var performance that relate to the use of multiple outer-loops, preconditioning, and the weak constraint are also explored. A systematic evaluation of the impact of the various components of the 4D-Var control vector (i.e. the initial conditions, surface forcing and open boundary conditions) is also presented. It is shown that correcting for uncertainties in the model initial conditions exerts the largest influence on the ability of the model to fit the available observations. Various important diagnostics of 4D-Var are also examined, including estimates of the posterior error, the information content of the observation array, and innovation-based consistency checks on the prior error assumptions. Using these diagnostic tools, we find that more than 90% of the observations assimilated into the model provide redundant information. This is a symptom of the large percentage of satellite data that are used and to some extent the nature of the data processing employed. This is the second in a series of three papers describing the ROMS 4D-Var systems.     

Picard iterations for a finite element shallow water equation model

ADCIRC, a finite element circulation model for shelves, coasts and estuaries, will be used for variational data assimilation. The nonlinear Euler–Lagrange (EL) problem will be solved using the iterated indirect representer algorithm. This algorithm makes such large, nonlinear but functionally smooth optimization problems feasible by iterating on linear approximations of the nonlinear problem (Picard iterations) and by making preconditioned searches in the “data subspace” at each iterate. Before solving the nonlinear EL using such Picard iterations, it essential that the iteration scheme be carefully examined within the framework of the nonassimilative or forward problem.The purpose of this paper is (1) to detail a Picard iteration procedure for ADCIRC, including the problematic bottom friction term; (2) to examine the ability of the iteration scheme to recover the nonlinear forward solution from deficient background fields; and (3) to present a study of different interpolation methods for reducing the memory/disk requirements of the iteration scheme. The iteration scheme is shown to be quite robust in its ability to recover the nonlinear solution from a variety of deficient background fields. A new cubic Hermitian interpolation method is shown to be a more effective alternative to standard linear interpolation for reducing memory/disk requirements, especially for high frequency overtides.     

Euler–Lagrange equations for the spectral element shallow water system

We present the derivation of the discrete Euler–Lagrange equations for an inverse spectral element ocean model based on the shallow water equations. We show that the discrete Euler–Lagrange equations can be obtained from the continuous Euler–Lagrange equations by using a correct combination of the weak and the strong forms of derivatives in the Galerkin integrals, and by changing the order with which elemental assembly and mass averaging are applied in the forward and in the adjoint systems. Our derivation can be extended to obtain an adjoint for any Galerkin finite element and spectral element system.We begin the derivations using a linear wave equation in one dimension. We then apply our technique to a two-dimensional shallow water ocean model and test it on a classic double-gyre problem. The spectral element forward and adjoint ocean models can be used in a variety of inverse applications, ranging from traditional data assimilation and parameter estimation, to the less traditional model sensitivity and stability analyses, and ensemble prediction. Here the Euler–Lagrange equations are solved by an indirect representer algorithm.     

A time‐split, forward‐backward numerical model for solving a nonhydrostatic and compressible system of equations

In this paper, we present a numerical procedure for solving a 2‐dimensional, compressible, and nonhydrostatic system of equations. A forward‐backward integration scheme is applied to treat high‐frequency and internal gravity waves explicitly. The numerical procedure is shown to be neutral in time as long as a Courant–Friedrichs–Lewy criterion is met. Compared to the leap‐frog‐scheme most models use, this method involves only two time steps, which requires less memory and is also free from unstable computational modes. Hence, a time‐filter is not needed. Advection and diffusion terms are calculated with a time step longer than sound‐wave related terms, so that extensive computer time can be saved. In addition, a new numerical procedure for the free‐slip bottom boundary condition is developed to avoid using inaccurate one‐sided finite difference of pressure in the surface horizontal momentum equation when the terrain effect is considered. We have demonstrated the accuracy and stability of this new model in both linear and nonlinear situations. In linear mountain wave simulations, the model results match the corresponding analytical solution very closely for all three cases presented in this paper. The analytical streamlines for uniform flow over a narrow mountain range were obtained through numerical integration of Queney's mathematical solution. It was found Queney's original diagram is not very accurate. The diagram had to be redrawn before it was used to verify our model results. For nonlinear tests, we simulated the famous 1972 Boulder windstorm and a bubble convection in an isentropic enviroment. Although there are no analytical solutions for the two nonlinear tests, the model results are shown to be very robust in terms of spatial resolution, lateral boundary conditions, and the use of the time-split scheme.     

Infrared atmospheric sounding interferometer data information content; instrument characterization and the impact of a priori information

The use of linear estimation for the study of the information content of a given satellite radiance data set for temperature and humidity profile retrievals is first reviewed. A particular formulation of the retrieval approach is then used to obtain an intrinsic characterisation of the Infrared Atmospheric Sounding Interferometer (IASI) data set, in terms of accuracy versus vertical resolution of retrieved profiles. The performance of the IASI instrument alone is analysed and compared to that of the currently‐used HIRS‐TOVS. The problem is then regularized by addition of a priori independent information to the initial data set. The potential use of IASI data for some particular choices of the a priori information associated with practical problems such as profile inversion or data assimilation for weather forecasting is analysed. The approach is finally used to derive an "empirical" objective framework to define the vertical discretization adapted to these problems.     

四维资料同化方法的特点分析和发展趋势

与把任意分布的观测值通过空间和时间插值分析到网格点的方法相比较 ,利用动力关系分析问题 ,其优势性十分明显。在一个能提供时间连续和动力耦合的模式预报方程中 ,有机地结合现在和过去的资料 ,即众所周知的四维资料同化。从目前经常应用的四维资料同化方法出发 ,详细分析了它们的特点和今后四维资料同化方法在预报模式中的应用前景。     

基于四维变分同化方法的南海中尺度涡后报实验

海洋中尺度涡在本质上是属于满足准地转平衡的大尺度运动,因此理论上,其在短时间内的运动将主要受到准地转平衡关系的约束,而外部强迫场的影响在短期内不会明显改变其运动特征。基于上述思想,我们提出了一种基于四维变分同化初始场的中尺度涡旋预报方案。为了检验该方案的可行性,本文使用区域海洋模式(regional ocean modeling system, ROMS)和其内建的增量强约束四维变分同化(incremental strong constraint four dimensional variational, I4D-Var)模块,建立了一个南海海洋同化模拟系统。首先,通过I4D-Var方法将AVISO卫星高度计资料同化到海洋数值模拟中,获得了理想的中尺度涡同化模拟结果。同化、模式模拟和观测三者的中尺度涡统计结果表明,该同化系统模拟的南海中尺度涡的路径、半径、海表高度异常和振幅等特征信息与AVISO(Archiving ValidationandInterpolationofSatelliteOceanographicData)观测结果高度吻合,同时在深度上的分析表明,涡旋对应的温度、盐度和密度均得到有效的调整。然后,将该同化系统的模拟结果做为初始场,对某一特定时段的南海中尺度涡进行了后报模拟和结果的定量化分析。通过比较后报模拟与观测资料中对应涡旋的海表面高度异常(sea surface height anomalies, SSHA)相关系数、涡心差距和半径绝对误差,证明该方案的中尺度涡后报时效至少可达10 d以上。后报实验结果验证了该中尺度涡预报方案的可行性,从而为中尺度涡的预报提供一定的理论基础和可行性方案。     

Ideal shocks in 2‐layer flow Part I: Under a rigid lid

Previous work on the classical problem of shocks in a 2‐layer density‐stratified fluid used either a parameterized momentum exchange or an assumed Bernoulli loss. We propose a new theory based on a set of viscous model equations. We define an idealized shock in two‐layer density stratified flow under a rigid lid as a jump or drop of the interface in which (1) the force balance remains nearly hydrostatic in the shock, (2) there is no exchange of momentum between the two layers except by pressure forces on the sloping interface, and (3) dissipative processes can be treated with a constant viscosity. We proceed in two steps. First, we derive a necessary condition for shock existence based on a requirement for wave steepening. Second, we formulate and solve a set of viscous model equations. Some results are the following: Shocks require strong layer asymmetry; one layer must be much faster and/or shallower than the other layer. The linearized equations describing the shock tails provide boundary conditions and a proof of shock uniqueness. It is possible to derive an analytical solution for weak shocks if the steepening condition is met. The weak shock solutions provide closed form expressions for the Bernoulli loss in each layer. Bernoulli losses are strongly concentrated in the expanding layer as the relative layer depth change is much larger in that layer. Bernoulli losses are independent of layer viscosity. A sudden cessation of shock existence is found for strong shocks when the possible end state migrates into the supercritical regime. Surprisingly, the new ideal shock theory compares well with a 2‐D, time‐dependent shallow water model (SWM) with a flux formulation, but with no viscous formulation. Both the Bernoulli drop and shock cessation condition agree quantitatively.     

基于地球系统模式的局地化粒子滤波器与集合卡尔曼滤波器同化实验

粒子滤波器(PF)是一种非常具有应用前景的非线性资料同化方法。但由于其算法本身存在的粒子退化问题,目前尚未被广泛地应用于大型地球物理模式。目前主流的集合同化系统仍然倾向于使用集合卡尔曼滤波器(EnKF)及其衍生方法。一种新近被提出的局地化粒子滤波器(LPF)在经典的粒子滤波器算法中引入局地化技术,可以使用较小的计算成本有效地避免退化问题,具有非常大的业务应用潜力。本文在全耦合的通用地球系统模式中开展了LPF和EnKF的同化实验,同化资料为模拟的卫星海表温度资料。着重考察了不同局地化参数对两种方法的不同影响,对比了局地化粒子滤波器与集合卡尔曼滤波器的同化效果差异。比较的结果表明,LPF的同化效果对于局地化参数的选择非常敏感,在使用最优局地化参数的条件下,LPF能达到与EnKF相当甚至优于后者的同化效果,并具有较大的改进空间。