Plasmonic metamaterials can exhibit a variety of physical optical properties that offer extraordinary nonlinear conversion efficiency for ultra-compact nanodevice applications. Furthermore, the optical-rectification effect from the plasmonic nonlinear metasurfaces (NLMSs) can be used as a compact source of deep-subwavelength thickness to radiate broadband terahertz (THz) signals. Meanwhile, a novel dual-mode metasurface consisting of a split-ring resonator (SRR) array and an epsilon-near-zero (ENZ) layer was presented to boost the THz conversion efficiency further. In this paper, to explore the mechanism of THz generation from plasmonic NLMSs, the Maxwell-hydrodynamic multiphysics model is adopted to investigate complex linear and intrinsic nonlinear dynamics in plasmonics. We solve the multiphysics model using the finite-difference time-domain (FDTD) method, and the numerical results demonstrate the physical mechanism of the THz generation processes which cannot be observed in our previous experiments directly. The proposed method reveals a new approach for developing new types of high-conversion-efficiency nonlinear nanodevices.