This study introduces the use of automatic differentiation (AD) to enable gradient-based optimization in two-dimensional, reactive flow simulations. A low-Mach formulation of the governing equations, expressed in primitive variables, is employed. The continuous Galerkin finite element method (CG-FEM), implemented in a semi-matrix-free manner, is used for spatial discretization. Operator splitting decouples the chemical and flow timescales, and a finite-rate chemistry model is utilized to represent the chemical kinetics. An implicit–explicit (IMEX) time integration scheme is applied, which facilitates stable AD over several milliseconds. The fully differentiable framework is realized in JAX, and verification is performed against an OpenFOAM (OF) solver. The validation cases include a one-dimensional, laminar, burner-stabilized flame and a two-dimensional, laminar, premixed slit flame. Two inverse modeling applications highlight the AD capabilities of the framework. In the first, a temperature boundary condition is inferred by minimizing discrepancies in flame shape based on two-dimensional heat release fields. In the second case, the parameters of a global two-step chemistry scheme are optimized to reproduce the flame tip position along the burner axis, as observed in experiments. The results confirm that the low-Mach, reactive CG-FEM solver delivers accurate forward simulations, while enabling stable AD-based inference and optimization.