Abstract: A multiscale derivative method is presented for solving a successive approximation linear quadratic regulator (SALQR) model for optimal in-situ bioremediation design. An efficient one-sided forward divided difference numerical derivatives calculation is implemented as the first stage of the method, which only requires assembling the right-hand-side vector of the linear systems of equations of the simulation model and performing backward substitution. The derivative calculation is reduced from O(N^3) to nearly O(N^2), where N is the number of non-Dirichlet state variables. A V-cycle multiscale derivatives approximation comprises the second stage, using coarser mesh derivatives to interpolate finer mesh derivatives. Implementing the numerical derivatives method in a case study with over 1600 state variables reduces computing time by more than two-thirds over the previous analytical derivatives method without loss of accuracy. The V-cycle multiscale derivatives approximation further reduces computing time by 29%, resulting in an overall 77% reduction compared to the previous analytical derivatives method. This reduction will be even greater for applications with more state variables, enabling solution of much larger-scale problems than previously possible.