Anion exchange membrane water electrolysis (AEMWE) is a promising route to low-cost green hydrogen, yet most existing studies remain at the single-cell scale, leaving the flow distribution, pressure management, and performance uniformity of industrial multi-cell stacks insufficiently understood. This study presents a modular design and optimization framework for multi-cell AEM electrolyzer stacks that combines a validated three-dimensional multiphysics computational fluid dynamics (CFD) model with the Taguchi method, Gaussian process regression (GPR) surrogate modeling, and an integrated sustainability assessment. A 24-cell stack model based on an Industrial Technology Research Institute (ITRI) prototype was developed in COMSOL Multiphysics and validated against experimental polarization data, yielding an average voltage deviation of 2.18%. Four geometric factors, namely the distribution plate manifold diameter, the electrolyte inlet/outlet diameter, the single-cell manifold diameter, and the single-cell inlet neck width, were investigated through an L9 orthogonal array to minimize pressure drop and maximize oxygen production. Analysis of variance was applied to quantify the relative contribution of each factor, and the optimized configuration was identified through signal-to-noise ratio analysis. A GPR surrogate trained on the L9 dataset was used to construct response surfaces and enable rapid prediction across the design space. An integrated assessment combining life cycle assessment (LCA), techno-economic analysis (TEA), and net carbon benefit (NCB) was performed to evaluate the global warming potential and levelized cost of hydrogen of the optimized stack, with Monte Carlo simulation of 10,000 iterations used to confirm the robustness of the estimates under input uncertainty. The proposed framework links geometric design, electrochemical performance, and sustainability metrics, providing practical guidance for the scale-up and commercialization of AEM electrolyzer stacks.