In this work, we present a transient LES approach using the Euler-Lagrange methodology and the Lattice-Boltzmann method, to successfully model and simulate interlinked bioreactor physics describing free surface hydrodynamics, multiphase mixing and mass transport in non-reacting and reacting systems. It is shown that the presence of reactions can result in a non-uniform spatially varying species concentration field, the magnitude and extent of which is directly related to the reaction rates and the underlying variations in the local volumetric mass transfer coefficient. Furthermore, we present a Gluconic acid optimization case study, which provides insights into the growth rates and process optimization over several hours of production. The computational approach presented in this work runs orders-of-magnitude faster than the conventional modeling tools and provides a unified framework for better understanding of free surface hydrodynamics, mixing, oxygen mass transfer and reaction kinetics to help with scale-up/scale-down decisions across a range of lengthscales.