Inverse problems governed by partial differential equations (PDEs) is a means to learn, from data, unknown or uncertain aspects of PDE-based mathematical models. PDE-based models often are constructed from scientific principles and are pervasive in science and engineering. The improvement of such a model can be accomplished by determining a ``best'' set of model parameters, or by reducing the uncertainty of the parameters that characterize the model. Ultimately, the goal is to improve the predictive capacity of such models and quantitatively understand the limitations of model-based predictions.In this dissertation, we describe algorithmic approaches for the Newton-based solution of large-scale computational inverse problems governed by PDEs. We present and motivate hierarchical matrix approximations of the Hessian, a key-component of Newton's method, as a means to exploit localized sensitivities of underlying elliptic PDE operators and which perform well for inverse problems with highly-informative data. To circumvent the computational challenges associated with generating hierarchical matrix approximations of matrix-free operators, such as the Hessian, we describe a local point spread function methodology, whose computational cost is independent of the problem discretization. It is numerically demonstrated that hierarchical matrix approximations of the Hessian can be effective for large-scale ice sheet inverse problems with highly informative data and thin ice sheets. Lastly, we present a scalable means to solve PDE- and bound-constrained optimization problems by an interior-point filter line-search strategy that leverages performant algebraic multigrid linear solvers. Bound constraints arise naturally in many inverse problems as a means to enforce sign-definiteness as dictated by e.g., physical principles. The inclusion of bound constraints does add particular computational challenges such as non-smooth complementarity conditions as part of the conditions for optimality. The Newton linear system solve strategy described in Chapter 4 is accelerated due to the inclusion of bound-constraints, wherein the performance of many other strategies is negatively impacted by such bound-constraints. A theoretical analysis is provided which demonstrates the algorithmic scaling of the approach. In addition, we demonstrate algorithmic and parallel scaling of the described approach by applying the framework to an example elliptic PDE- and bound-constrained problem