As modern field effect transistors (FETs) scale down into the sub-10-nm regime, the approaching end of CMOS scaling becomes evident. Recently, we have performed a series of three-dimensional, fully quantum-mechanical, charge-self-consistent transport simulations and optimizations to show [1] that all FETs, irrespective of their channel material, will reach a fundamental downscaling limit around 4-5 nm technology node, due to unacceptably high levels of thermally induced errors. To investigate alternatives to FET/CMOS scaling, many emerging devices and technologies have been under active research in recent years. Among them, oxide-based memristors have become one of the strongest candidates to replace flash, and possibly DRAM and SRAM in the near future [2]. Moreover, memristors have a high potential to enable beyond-CMOS technology advances in novel architectures for high performance computing [3]. However, despite the spectacular progress in experimental demonstration of oxide memristors, their microscopic transport theory and even some key operational principles, essential to advance experimental progress, remain uncertain [4]. To facilitate the understanding of physical switching mechanisms and accelerate experimental development in oxide memristors, we have developed a 3D fully-coupled electrical and thermal transport model [5] that solves simultaneously the time-dependent continuity equations for electrons, holes, and oxygen vacancies, the time-dependent heat equation including Joule heating sources, and the Poisson equation for all the charged species. The model captures all the most important processes that drive memristive switching, including field drift due to the electrical potential gradient, Fick diffusion due to the concentration gradient, and the Soret effect due to the temperature gradient. To properly simulate the movement of oxygen vacancies, the rigid point ion model [5] is used, in which the vacancy velocity exponentially depends on the temperature and the electric field when the field exceeds a critical value. The fully-coupled equations and relevant physical models are implemented in CHARON, a multi-dimensional parallel TCAD code developed at Sandia National Labs. CHARONis built upon the open-source Trilinos packages [6], which support finite element and finite volume discretization methods, state-of-the-art nonlinear and linear solvers, MPI parallel capability, and many more. It can simulate arbitrary device geometry and supports multi-physics capability (i.e., allowing for solving different equations in different regions). The coupled model has been applied to simulate the switching process in a 3D filamentary tantalum oxide memristor shown in Fig. 1(a). The layer material and thickness come from an experimental device published in Ref. [7]. Making use of the multi-physics capability in CHARON, we solve the coupled Poisson, electron, oxygen vacancy, and heat equations in the active Ta and conduction filament (CF) regions, with vacancies confined in Ta and CF regions, while in the Pt and Ta2O5 regions, due to absence of mobile vacancies, the coupled Poisson, electron, and heat equations are solved, with Pt treated as a highly doped n-type semiconductor using immobile dopants, and Ta2O5 as an undoped semiconductor. Extensive simulations show that both field drift and thermal processes play crucial roles in determining the switching dynamics of the bipolar TaOx device. Specifically, (i) during OFF switching, when a sufficient negative voltage is applied to the top Pt electrode, Joule heating causes the increase of temperature with the hottest spot located in the CF, and consequently, oxygen vacancies move away from CF into Ta under the influence of both heating and field, which leads to a vacancy density gap and possible depletion of vacancies in the entire CF at high negative voltage, resulting in the OFF state; (ii) during ON switching, a sufficient positive voltage causes Joule heating and the increase of temperature, and increased temperature causes vacancies moving from Ta into CF to refill the vacancy gap, resulting in the ON state. Interestingly, we show that the thermally activated field-dominated switching in the considered devices can lead to more than three orders of magnitude speed-up in switching time, with only 50% voltage/current increase. Simulated resistances during OFF and ON switching as a function of time and voltage show very good qualitative agreement with experimental measurement in Ref. [7]. Furthermore, under a triangular voltage sweep, the simulated current-voltage hysteresis curve is shown in Fig. 1(b), which also shows good qualitative agreement with the measured data in Ref. [8]. The generality of the proposed model allows to simulate and investigate the underlying physical switching mechanisms in a wide spectrum of oxide memristors ranging from field-dominant, field-thermal-driven to thermal-dominant memristive devices. Figure 1