In this paper, a high-accuracy H 1 -Galerkin mixed finite element method (M-FEM) for strongly damped wave equation is studied by linear triangular finite element. By constructing a suitable extrapolation scheme, the convergence rates can be improved from O(h) to O(h 3 ) both for the original variable u in H 1 (Ω) norm and for the actual stress variable p = ∇u t in H(div;Ω) norm, respectively. Finally, numerical results are presented to confirm the validity of the theoretical analysis and excellent performance of the proposed method. 611 As we know, H 1 -Galerkin MFEM was firstly proposed in [20] for solving parabolic problem. It has the same convergence rates as the standard MFEMs, but without the requirement of LBB consistency condition and quasi-uniformity requirement on the subdivision meshes, which allow the approximation spaces to be chosen freely. Recently, there are some extensive studies on H 1 -Galerkin MFEMs for parabolic partial integrodifferential equation [21], second order hyperbolic equation [22], nonlinear parabolic equation in porous medium flow [2], pressure equation in compressible porous medium flow [3], pseudo-hyperbolic equation [18], and Sobolev equation [26], etc. As for the strongly damped wave equation (1.1), the optimal order error estimates of H 1 -Galerkin MFEM were investigated in [17]. Recently, though there are many excellent work on superconvergence analysis of MFEMs [4, 5, 13, 16], there are not much work devoted to the development of superconvergence analysis of H 1 -Galerkin MFEM. Via integral equalities techniques and a postprocessing technique, in [32] and [31], Tripathy and Sinha derived superconvergence estimates of H 1 -Galerkin MFEM by using rectangular Raviart-Thomas finite element for elliptic equation and parabolic problem, respectively. Later, [28,30] studied superconvergence of H 1 -Galerkin MFEM by the linear triangular finite element and nonconforming quasi-Wilson finite element, respectively, and improved the approximation accuracy from O(h) to O(h 2 ) both for u in H 1 (Ω) norm and the actual stress variable p = ∇u t in H(div;Ω) norm.On the other hand, Richardson extrapolation is often used to increase the convergence rate of numerical solution by combining the two numerical solutions obtained on different meshes. Practical application of this approach to FEM can be found in [8,11,12]. As to MFEM, [6,9,15] gave the asymptotic error expansions of the lowest rectangular Raviart-Thomas finite element so as to improve the approximation accuracy in L 2 (Ω) norm for second order elliptic problem, parabolic type integro-differential equation and eigenvalue problem with integral identities and Richardson extrapolation techniques. Moreover, [33] studied the extrapolation of the Nédélec finite element for Maxwell equation.[35] also investigated the asymptotic error expansions for the Stokes eigenvalue problem by Bernadi-Raugel finite element.However, to our best knowledge, there are few reports on the study of extrapolations of H 1 -Galerkin MFEM except for our early wo...