Gold nanohole arrays (GNAs) are hybrid metal/dielectric metasurfaces consisting of periodically arranged cylindrical air holes drilled in a gold film. Thanks to their morphology, GNAs can support both localized and propagating surface plasmons, providing high suitability for biosensing applications based on surface plasmon resonance (SPR). However, both optical properties and SPR performances are significantly influenced by the GNA geometrical parameters and, as a consequece, in order to obtain a valuable and controlled response a precise tailoring of these features is needed. For this purpose, we investigated the impact of the array's pitch, nanohole radius, and gold thickness, while also considering nanofabrication restrictions. The study was made on a square GNA employing three complementary computational approaches. A design routine based on a custom particle swarm optimization (PSO) algorithm implemented in Ansys Lumerical FDTD was used to obtain, for a fixed gold thickness (spanning from 40 nm to 100 nm), the array’s pitch, and nanohole radius corresponding to a GNA structure characterized by an optimized and precisely tuned optical response. Specifically, by mimicking a real SPR measurement, the algorithm provides the geometrical parameters that maximize the coupling of the source with the main localised plasmonic mode. In addition, considering SPR instrumental constraints, the tuning of the reflectance minima at a wavelength of 770 nm was imposed. Starting from the structures given by the PSO, a consequent custom sweep, also implemented in Ansys Lumerical FDTD, was used to refine the optical response in terms of the spectral position for the reflectance minimum as a function of the GNA geometrical parameters. For the spanned gold thicknesses, the reflectance minima positions exhibit a cubic-like dependence on both the array pitch and nanohole radius, consistently with the Bragg condition. To support and better explain these effects in terms of scattering properties of cylindrical voids in a metallic film with finite thickness, a scattering matrix formalism was applied to the analysis. It also enables a comparison of the electric field profiles with the field expansion given by the FDTD. In addition, to understand and link the results to the excitation and plasmonic modes dispersion, a band structure analysis was carried out, again using Ansys Lumerical FDTD, providing a more photonic crystal-like model of the GNA metasurface.