AbstractIncremental 4D‐Var is a data assimilation algorithm used routinely at operational numerical weather prediction (NWP) centres worldwide. The algorithm solves a series of quadratic minimization problems (inner‐loops) obtained from linear approximations of the forward model around nonlinear trajectories (outer‐loops). Since most of the computational burden is associated with the inner‐loops, many studies have focused on developing computationally efficient algorithms to solve the least‐square quadratic minimization problem, in particular through time parallelization. This paper presents the first implementation and testing of a recently proposed method for parallelizing incremental 4D‐Var, the Randomized Incremental Optimal Technique (RIOT), which replaces the traditional sequential conjugate gradient (CG) iterations in the inner‐loop of the minimization with fully parallel randomized singular value decomposition (RSVD) of the preconditioned Hessian of the cost function. RIOT is tested using the standard Lorenz‐96 model (L‐96) as well as two realistic high‐dimensional atmospheric source inversion problems based on aircraft observations of black carbon concentrations. A new outer‐loop preconditioning technique tailored to RSVD was introduced to improve convergence stability and performance. Results obtained with the L‐96 system show that the performance improvement from RIOT compared to standard CG algorithms increases significantly with nonlinearities. Overall, in the realistic black carbon source inversion experiments, RIOT reduces the wall‐clock time of the 4D‐Var minimization by a factor of 2 to 3, at the cost of a factor of 4 to 10 increase in energy cost due to the large number of parallel cores used. Furthermore, RIOT enables reduction of the wall‐clock time computation of the analysis‐error covariance matrix by a factor of 40 compared to a standard iterative Lanczos approach. Finally, as evidenced in this study, implementation of RIOT in an operational NWP system will require a better understanding of its convergence properties as a function of the Hessian characteristics and, in particular, the degree of freedom for signal of the inverse problem.