Full-waveform inversion (FWI) has become a powerful tool in inverting subsurface geophysical properties. The estimation of uncertainty in the resulting Earth models and parameter trade-offs, although equally important to the inversion result, has however often been neglected or became prohibitive for large-scale inverse problems. Theoretically, the uncertainty estimation is linked to the inverse Hessian (or posterior covariance matrix), which for massive inverse problems becomes impossible to store and compute. In this study, we investigate the application of the square-root variable metric (SRVM) method, a quasi-Newton optimisation algorithm, to FWI in a vector version. This approach allows us to reconstruct the final inverse Hessian at an affordable storage memory cost. We conduct SRVM based elastic FWI on several elastic models in regular, free-surface and practical cases. Comparing the results with those obtained by the state-of-the-art L-BFGS algorithm, we find that the proposed SRVM method performs on a similar, highly-efficient level as L-BFGS, with the advantage of providing additional information such as the inverse Hessian needed for uncertainty quantification.