Uncertainty quantification (UQ) in simulation-based science seeks to estimate magnitudes of errors in the model inputs (prior uncertainties), assess - through formal inversion - to which extent observational constraints, along with their instrument and representation errors, can reduce these errors (posterior uncertainties), and propagate the (prior or posterior) error covariances of model inputs to key target quantities of interest (climate indices). Uncertain inputs in Earth system models (EaSMs) span a very high-dimensional space (e.g., 10^8), consisting of 3D initial condition or model parameter fields (capturing model error), as well as 2D surface or basal boundary conditions, sometimes with significant time-variation (forcing uncertainties). Derivative-based methods are a powerful tool to make otherwise computationally prohibitive UQ inferences tractable. Adjoint models are a highly efficient method for providing gradient information for optimization or sensitivity information for science, while Hessian-based methods provide measures of posterior uncertainty covariances. Obtaining derivative codes for complex EaSMs remains a challenge, but which has been successfully tackled through combined use of algorithmic differentiation (e.g., using OpenAD) and solution of the adjoint equations via efficient solver packages. Here we provide three example from ice sheet, ocean, and biogeochemical modeling.