The coupling between plate tectonics and mantle convection is generally investigated with numerical resolution methods that prescribe the number of plates a priori. We have designed a new MACMA model (MultiAgent Convecting MAntle) that simulates time-dependent plate tectonics in a 2D cylindrical geometry with evolutive plate boundaries. Plate velocities are computed using a local force balance and explicit parameterizations are used to treat plate boundary processes such as trench migration, subduction initiation, continental breakup and plate suturing. The system's geometry and thermal state are subsequently updated at all times. In our model, the number of plates is not imposed but emerges naturally. This approach is based on multiagent systems in thermal and mechanical interaction. Our model involves four types of agents (convection cells, lithospheric plates, continents and plate boundaries) that collect information from their environment in order to make decisions controlled by behavior laws. Our approach has two goals: (1) to test how empirically- and analytically-determined rules for a specific agent's behavior affect plate dynamics as a whole, and (2) to investigate how mantle temperature evolution is influenced by evolving surface plate tectonics. We obtain predictions for driving forces and corresponding plate velocities that are in good agreement with observations. Our model highlights the key role played by the sea floor age distribution, and thus by local structure changes (e.g., creation of a new spreading ridge), on the short-term heat flux evolution. Moreover, the present-day thermal and geometrical state of the planet corresponds to a decreasing heat flow, as suggested by recent seafloor age reconstructions. In the long term, Earth's thermal history is classically studied using scaling laws that link surface heat loss to the temperature and viscosity of the convecting mantle. When such a parameterization is used in the global heat budget of the Earth to integrate the mantle temperature backwards in time, a runaway increase of temperature is obtained, leading to the so-called "thermal catastrophe". When investigating Earth's cooling rate using explicit tectonic mechanisms, the very low computational cost of our model allows us to study the effect of a wide range of input parameters on the long-term thermal evolution of the system. For Earth-like parameters, an average cooling rate of 60~K per billion years is obtained, which is consistent with petrological and rheological constraints. Two time scales arise in the evolution of the heat flux: a linear long-term decrease and high-amplitude short-term fluctuations due to tectonic rearrangements. We show that the viscosity of the mantle is not a key parameter in the thermal evolution of the system. The cooling rate of the Earth depends mainly on its ability to replace old insulating seafloor by young thin oceanic lithosphere. Therefore, the main controlling factors are parameters such as the resistance of continental lithosphere to breakup or the critical age for subduction initiation. We infer that simple convective considerations alone cannot account for the complex nature of mantle heat loss and that tectonic processes dictate the thermal evolution of the Earth.
