The In-Situ Upgrading (ISU) of bitumen and oil shale is a very challenging process to model numerically because a large number of components need to be modelled using a system of equations that are both highly non-linear and strongly coupled. In addition to the transport of heat by conduction and convection, and the change of properties with varying pressure and temperature, these processes involve transport of mass by convection, evaporation, condensation and pyrolysis chemical reactions. The behaviours of these systems are difficult to predict as relatively small changes in the material composition can significantly change the thermophysical properties. Accurate prediction is further complicated by the fact that many of the inputs needed to describe these processes are uncertain, e.g. the reaction constants and the temperature dependence of the material properties. The large number of components and chemical reactions involves a non-linear system that is often too large for full field simulation using the Fully Implicit Method (FIM). Operator splitting (OS) methods are one way of potentially improving computational performance. Each numerical operator in a process is modelled separately, allowing the best solution method to be used for the given numerical operator. A significant drawback to the approach is that decoupling the governing equations introduces an additional source of numerical error, known as splitting error. Obviously the best splitting method for modelling a given process is the one that minimises the splitting error whilst improving computational performance over that obtained from using a fully implicit approach. Although operator splitting has been widely used for the modelling of reactive-transport problems, it has not yet been applied to models that involve the coupling of mass transport, heat transfer and chemical reactions. One reason is that it is not clear which operator splitting technique to use. Numerous such techniques are described in the literature and each leads to a different splitting error, which depends significantly on the relative importance of the mechanisms involved in the system. While this error has been extensively analysed for linear operators for a wide range of methods, the results observed cannot be extended to general non-linear systems. It is therefore not clear which of these techniques is most appropriate for the modelling of ISU. Analysis using dimensionless numbers can provide a useful insight into the relative importance of different parameters and processes. Scaling reduces the number of parameters in the problem statement and quantifies the relative importance of the various dimensional parameters such as permeability, thermal conduction and reaction constants. Combined with Design of Experiments (DOE), which allows quantification of the impact of the parameters with a minimal number of numerical experiments, dimensionless analysis enables experimental programmes to be focused on acquiring the relevant data with the appropriate accuracy by ranking the different parameters controlling the process. It can also help us design a better splitting method by identifying the couplings that need to be conserved and the ones that can be relaxed. This work has three main objectives: (1) to quantify the main interactions between the heat conduction, the heat and mass convection and the chemical reactions, (2) to identify the primary parameters for the efficiency of the process and (3) to design a numerical method that reduces the CPU time of the simulations with limited loss in accuracy. We first consider a simplified model of the ISU process in which a solid reactant decomposes into non-reactive gas. This model allows us to draw a parallel between the in-situ conversion of kerogen and the thermal decomposition of polymer composite when used as heat-shield. The model is later extended to include a liquid phase and several reactions. We demonstrate that a ISU model with nf fluid components, ns solid components and k chemical reactions depends on 9+k*(3+nf+ns-2)+8nf+2ns dimensionless numbers. The sensitivity analysis shows that (1) the heat conduction is the primary operator controlling the time scale of the process and (2) the chemical reactions control the efficiency of the process through the extended Damköhler numbers, which quantify the ratio of chemical rate to heat conduction rate at the heater temperature for each reaction in the model. In the absence of heat loss and gravity effects, we show that the ISU process is most efficient at a heater temperature for which the minimum of the extended Damköhler numbers of all reactions included in the model was between 10 and 20. For the numerical method, the standard Iterative Split Operator (ISO) does not perform well due to many convergence failures, whereas the standard Sequential Split Operator (SSO) and the Strang-Marchuk Split Operator (SMSO) give large discretization errors. We develop a new method, called SSO-CKA, which has smaller discretization error. This method simply applies SSO with three decoupled operators: the heat conduction (operator $C$), the chemical reactions (operator $K$) and the heat and mass convection (operator $A$), applied in this order. When we apply SSO-CKA with the second-order trapezoidal rule (TR) for solving the chemical reaction operator, we obtain a method which generally gives smaller discretization errors than FIM. We design an algorithm, called SSO-CKA-TR-AIM, which is faster and generally more accurate than FIM for simulations with a kinetic model including a large number of components that could be regrouped into a small number of chemical classes for the advection and heat conduction operator. SSO-CKA works best for ISU models with small reaction enthalpies and no other reaction than pyrolysis reactions, but can give a large discretization method for ISU models with non-equilibrium reactions.