Solving the time-dependent quantum many-body Schrödinger equation is a challenging task, especially for states at a finite temperature, where the environment affects the dynamics. Most existing approximating methods are designed to represent static thermal density matrices, 1D systems, and/or zero-temperature states. In this work, we propose a method to study the real-time dynamics of thermal states in two dimensions, based on thermofield dynamics, variational Monte Carlo, and neural-network quantum states. To this aim, we introduce two novel tools: (i) a procedure to accurately simulate the cooling down of arbitrary quantum variational states from infinite temperature, and (ii) a generic thermal (autoregressive) recurrent neural-network (ARNNO) Ansatz that allows for direct sampling from the density matrix using thermofield basis rotations. We apply our technique to the transverse-field Ising model subject to an additional longitudinal field and demonstrate that the time-dependent observables, including correlation operators, can be accurately reproduced for a 4×4 spin lattice. We provide predictions of the real-time dynamics on a 6×6 lattice that lies outside the reach of exact simulations.