Abstract Equation discovery methods, such as symbolic regression, show great promise to generate parameterizations of biogeochemical processes in an objective data‐driven manner, yet remain untested in ocean biogeochemistry. Here, we apply symbolic regression to a state‐of‐the‐art ocean biogeochemical model, using it as a surrogate data set to rediscover an empirical equation used to calculate colloidal iron in the model. We introduce a robustness metric combining R 2 (global pattern reproductio