Finite Element (FE) modeling often requires unavoidable simplifications or assumptions due to a lack of experimental data, modeling complexity, or non-affordable computational cost. One such simplification is modeling corrosion phenomena or material properties, which are usually assumed to be uniform throughout the structure. However, e.g., corrosion has a local nature and severe consequences on the behavior of steel structures that should not be overlooked. To improve the current numerical modeling techniques in aging steel bridges, this paper proposes a Gaussian Copula-based Bayesian Network (GCBN) approach to model the spatial variability of structural element properties. Accordingly, a study of the automatic Bayesian network generation process is first conducted. Subsequently, the methodology is applied to a severely damaged riveted steel bridge built in 1897. The results show that the methodology has excellent flexibility for generating properties variability in FE models at a low computational cost, thus ensuring its practical feasibility and robustness for accurate numerical modeling.