A large number of new measurements with the activation technique were performed for (n,2n) and neutron-induced Z = 1,2 reaction cross sections on the stable molybdenum isotopes in the energy range from 13.5 to 21 MeV. First results were obtained for the 92 Mo(n,2n) Mo(n,x) 97 Nb m , and 100 Mo(n,α) 97 Zr reactions, above 16 MeV. A significant number of high-accuracy 14 MeV measurements were performed which are in good agreement with the measurements above 16 MeV for reactions studied in both energy ranges. The rather complete database for the molybdenum isotopes was analyzed with two different sets of consistent model calculations: a local and a global approach. The global approach (a blind calculation with the TALYS code) provides a good overall description of the dominant reaction channels, although the (n,α) reactions for the heavy isotopes are overpredicted. The local approach (an adjusted calculation with the STAPRE-H code) describes the shapes and magnitudes of the excitation functions well from the reaction thresholds up to 21 MeV using a consistent parameter set, which was optimized based on all experimental information for the nuclei at hand and their immediate neighbors. The agreement between experimental and calculated data is, in general, good both at the maxima and at the tails of the excitation functions, and both for total activation cross sections of a particular channel and for cross sections leading to isomers, showing the viability of the level densities, the optical models, and the γ widths. Comparison of the two model calculations with the data indicates the relevance of an appropriate treatment for preequilibrium (PE) α-particle emission for the description of the data above 14 MeV. Comparison between the model calculations shows largely different PE deuteron emission contributions to the total ( Z = 1, A = 1) cross sections with an additional marked difference in energy dependence. This suggests that emission spectra around 20 MeV are required to establish the magnitude of the PE deuteron emission contribution to this process. New γ -ray strength functions were established by verification against average (n,γ ) data and were demonstrated to give good agreement with the measured isomer production cross sections.