Although in FWI the gradient optimisation-based methods have been well developed, the uncertainty estimation methods are also essential but still left behind. The evaluation of FWI uncertainty involves the Hessian-related posterior covariance matrix, which is prohibitive to compute and store for practical problems. To confront this issue, we propose to run FWI with square-root-variable-metric (SRVM), a quasi-Newton method similar to L-BFGS, in a matrix-free vector version. The vector-version SRVM is memory-affordable even for large-scale problems. The size of one SRVM vector is identical to that of the parameter model, and the number of SRVM vectors equals the total iteration number. We validate SRVM within acoustic FWI, with L-BFGS for reference. After the SRVM-based FWI converges, we can access the approximated inverse data-misfit Hessian and then the posterior covariance from the stored SRVM scalar and vector series. To reconstruct the inverse Hessian more efficiently, we factorise its eigenvalues and eigenvectors via randomised singular value decomposition (RSVD). We can qualify the inversion uncertainty through the posterior standard deviation, the 2D prior/posterior random samplings. We demonstrate our methods with numerical examples.