Abstract: The late-time component of high-altitude electromagnetic pulse (HEMP E3) is a serious electromagnetic threat to modern power grids because it can induce quasi-DC geomagnetically induced current (GIC) in long transmission lines. When these currents intrude into power transformers through grounded neutral points, they cause DC magnetic bias and half-wave saturation of the iron core. This saturation leads to excitation-current distortion, increased electromagnetic losses, abnormal stray losses, and possible local overheating. Accurate evaluation of the electromagnetic and thermal responses of transformers under HEMP E3-induced GIC conditions is therefore important for assessing transformer operating safety and power-grid resilience. This study aims to develop a two-dimensional electromagnetic-fluid-thermal coupling calculation framework for analyzing the electromagnetic response, loss distribution, oil natural-convection cooling process, and transient temperature-rise behavior of oil-immersed transformers under late-time HEMP excitation. To achieve this objective, a field-circuit coupled electromagnetic model was established for a three-phase transformer with a grounded neutral point. Under the low-frequency magneto-quasi-static assumption, the edge-based finite element method was adopted to solve the magnetic vector potential equation. The winding circuit equations were coupled with the electromagnetic field equations through the winding current density and the induced electromotive force. The nonlinear magnetization behavior of the iron core was described by the B-H magnetization curve of silicon steel, and a relaxation strategy was introduced during the iterative update of current density to improve numerical stability under DC magnetic bias. After the electromagnetic solution was obtained, the hysteresis loss, eddy-current loss, winding Joule loss, and stray loss were calculated and mapped into the thermal model as non-uniform volumetric heat sources. For the fluid-thermal calculation, the transformer oil was treated as an incompressible Newtonian fluid under natural oil-circulation cooling conditions. The Boussinesq approximation was used to describe buoyancy-driven flow caused by temperature-dependent density variation. The Brinkman penalization model was introduced to distinguish the fluid region from the solid regions in the unified computational domain. The Navier-Stokes equations and energy conservation equation were solved in the oil region, while the heat conduction equation was solved in the iron core, windings, and other solid regions. A segregated transient solution procedure was used to update the velocity, pressure, and temperature fields at each time step. Transient simulations were then carried out on a two-dimensional three-phase oil-immersed transformer model, and the calculated results were compared with those obtained from COMSOL. The results show that, under the normal operating condition without GIC injection, the calculated primary and secondary voltages remain sinusoidal and maintain the expected phase relationship and voltage ratio. The primary current exhibits typical excitation-current distortion caused by nonlinear magnetic behavior, while the positive and negative half cycles remain symmetric, indicating the absence of DC magnetic bias. The temperature-rise curve at the selected monitoring point agrees well with the COMSOL result. Under the DC magnetic bias condition with GIC injection, the calculated voltage waveforms are also consistent with the COMSOL results in phase, amplitude, and waveform shape. The proposed model can capture the dynamic variation of the main magnetic flux under nonlinear saturation and can describe the heat accumulation caused by electromagnetic losses and heat dissipation through transformer-oil natural convection. The results demonstrate that the proposed two-dimensional electromagnetic-fluid-thermal coupling method has good convergence and computational accuracy for analyzing transformer responses under HEMP E3-induced GIC conditions. The method provides a useful numerical tool for evaluating electromagnetic loss distribution and transient temperature-rise behavior of oil-immersed transformers under natural oil-circulation cooling conditions. The present model is mainly intended for coupling-algorithm verification and temperature-rise trend analysis. For engineering-level prediction of local hot spots in full three-dimensional transformers, further extension to three-dimensional structures, forced-oil or directed-oil cooling boundaries, cooling-channel details, and more refined thermal boundary conditions is still required.