Summary: Inverse problems in ocean acoustics are generally solved by means of matched-field inversion in combination with meta-heuristic global search algorithms. Recently adjoint-based optimal nonlocal boundary control has been proposed as a possible complement or alternative to the traditional geoacoustic inversion methods. Especially for the assessment of a shallow water environment where field data are recorded on non-fully populated hydrophone arrays there is an inherent need for observation over a broad range of frequencies, as has been widely experimentally demonstrated. This article presents an analytic multiple frequency adjoint formulation for estimating an acoustically equivalent model of the seabed based on local impedance boundary conditions. It is intended as a benchmark solution for an extension of the approach to more sophisticated boundary conditions by means of automatic differentiation. The approach, coming from control theory, aims at inverting the impedance boundary condition along the water-sediment interface via measurements of the acoustic propagation in the water column. The impedance boundary condition thus plays the role of the control parameter of the complex pressure field in the waveguide. The developed multiple frequency adjoint method of optimal control provides an explicit representation of the gradient of the matched-field oriented cost function, and the article further proposes an analytical second-order adjoint formulation which permits the computation of the product of the Hessian matrix by a vector. By means of example inversion results the feasibility of the multi-frequency approach is discussed and the effect of a limited set of receiving locations is illustrated.