We investigate scalable multilevel domain decomposition algorithms for solving inverse elliptic problems formulated as optimization problems constrained by partial differential equations. To solve these optimization problems, we employ a fully coupled Lagrange-Newton-Krylov-Schwarz algorithm. One of the key steps in the algorithm is the Jacobian preconditioning, for which we study and compare four types of two-level domain decomposition methods. Our numerical results show that the algorithms work well for different types of observations in terms of the accuracy of the solution, and some of the algorithms scale better than the others when the number of processors is large. We also study and report the sensitivity of the algorithms with respect to the jumps of the coefficients, the level of noise in the observed data, the size of the computational domain, the size of the mesh, and the number of processors.