Bubble reactors using molten metal alloys (e.g, nickel-bismuth and copper-bismuth) with strong catalytic activity for methane decomposition are an emerging technology to lower the carbon intensity of hydrogen production. Methane decomposition occurs non-catalytically inside the bubbles and catalytically at the gas-liquid interface. The reactor performance is therefore affected by the hydrodynamics of bubble flow in molten metal, which determines the evolution of the bubble size distribution and of the gas holdup along the reactor height. A reactor model is first developed to rigorously account for the coupling of hydrodynamics with catalytic and non-catalytic reaction kinetics. The model is then validated with previously reported experimental data on methane decomposition at several temperatures in bubble columns containing a molten nickel-bismuth alloy. Next, the model is applied to optimize the design of multitubular catalytic bubble reactors at industrial scales. This involves minimizing the total liquid metal volume for various tube diameters, melt temperatures, and percent methane conversions at a specified hydrogen production rate. For example, an optimized reactor consisting of 891 tubes, each measuring 0.10 m in diameter and 2.11 m in height, filled with molten Ni0·27Bi0.73 at 1050 °C and fed with pure methane at 17.8 bar, may produce 10,000 Nm3.h−1 of hydrogen with a methane conversion of 80% and a pressure drop of 1.6 bar. The tubes could be heated in a fired heater by burning either a fraction of the produced hydrogen, which would prevent CO2 generation, or other less expensive fuels.