In this study we have employed statistical methods to efficiently design experiments and analyze output of an ocean general circulation model that uses an isopycnal mixing parameterization. Full ranges of seven inputs are explored using 51 numerical experiments. Fifteen of the cases fail to reach satisfactory equilibria. These are attributable to numerical limitations specific to the isopycnal model. Statistical approximating functions are evaluated using the remaining cases to determine the dependency of each of the six scalar outputs on the inputs. With the exception of one output, the approximating functions perform well. Known sensitivities, particularly the importance of diapycnal (vertical) eddy diffusivity and wind stress, are reproduced. The sensitivities of the model to two numerical constraints specific to the isopycnal parameterization, maximum allowable isopycnal slope and horizontal background eddy diffusivity, are explored. Isopycnal modelling issues, convection reduction and the Veronis effect, are examined and found to depend crucially on the isopycnal modelling constraints.