Unsteady aerodynamic phenomena such as gusts and flutter impose significant loads on aircraft and must be predicted accurately during design, certification, and optimisation. Although higher-fidelity methods like solving the unsteady Reynolds-Averaged Navier-Stokes (URANS) equations or large-eddy simulations can capture complex unsteady behaviour, including dynamic non-linearities, their computational cost is prohibitive for many-query applications. Machine learning (ML)–based surrogate models offer a promising alternative, combining the ability to capture non-linear relationships with fast evaluation once trained. For computational fluid dynamics the spatial domain is discretised as a mesh. Graph neural networks (GNNs) are a ML method that can exploit the mesh structure as an inductive bias [1, 2]. This work introduces an URANS-based dataset containing multiple trajectories of unsteady pitch and heave motions of a two-dimensional airfoil with varying amplitudes, frequencies, and initial conditions. The goal is to find a surrogate model that predicts surface quantities (pressure and skin friction) given a steady initial flow solution and an exogenous input signal that describes the airfoil’s motion in time. Initial investigations that extend the approach from [2] with exogenous inputs identified the error accumulation of the autoregressive model as the main limitation for long roll-outs in time. Several approaches to mitigate this problem are investigated. Firstly, architectural adaptations such as the integration of temporal ML components (Gated Recurrent Units or temporal attention) inside the GNN structure are explored. Secondly, training methodologies are evaluated. Finally, the approach is reformulated as a neural ordinary differential equation (Neural ODE) [3], combining strengths of classical numerical methods with ML. While the first two approaches lead to some improvements, preliminary results on the latter approach look most promising.