Ensemble Kalman methods were initially developed to solve nonlinear data assimilation problems in oceanography, but are now popular in applications far beyond their original use cases. Of particular interest is climate model calibration. As hybrid physics and machine-learning models evolve, the number of parameters and complexity of parameterizations in climate models will continue to grow. Thus, robust calibration of these parameters plays an increasingly important role. We focus on learning climate model parameters from minimizing the misfit between modeled and observed climate statistics in an idealized setting. Ensemble Kalman methods are a natural choice for this problem because they are derivative-free, scalable to high dimensions, and robust to noise caused by statistical observations. Given the many variants of ensemble methods proposed, an important question is: Which ensemble Kalman method should be used for climate model calibration? To answer this question, we perform systematic numerical experiments to explore the relative computational efficiencies of several ensemble Kalman methods. The numerical experiments involve statistical observations of Lorenz-type models of increasing complexity, frequently used to represent simplified atmospheric systems, and some feature neural network parameterizations. For each test problem, several ensemble Kalman methods and a derivative-based method "race" to reach a specified accuracy, and we measure the computational cost required to achieve the desired accuracy. We investigate how prior information and the parameter or data dimensions play a role in choosing the ensemble method variant. The derivative-based method consistently fails to complete the race because it does not adaptively handle the noisy loss landscape.