In contrast to many other types of cancer, the incidence of liver cancer, specifically hepatocellular carcinoma (HCC), is on the rise. For most patients, surgical intervention is not a viable option, leaving them reliant on chemotherapy treatments, particularly transarterial chemoembolization (TACE), for relief. Our study aims to understand how these treatments function within the liver and their impact on tumor growth. Building upon existing research, we model the flow and transport of chemotherapy drugs and embolic agents in the liver using the miscible displacement equations. Utilizing CT images from liver cancer patients, we extract a 1D centerline of the hepatic vascular structures that deliver blood to the tumors, and then construct a 3D mesh from the liver segmentations. We employ the singularity subtraction technique to create a finite element model for the flow of blood in the liver, specifically focusing on areas affected by the TACE treatment. We extend the singularity subtraction technique to the time-dependent advection-diffusion equation to model the concentration of chemotherapy drugs in the liver and tumors. We first solve the time-dependent non-conservative advection-diffusion equation using the finite element method. To address instabilities arising when the model is advection dominated, we then utilize the discontinuous Galerkin method to solve the time-dependent conservative advection-diffusion equation. We couple the models for blood flow following the injection of an embolic agent with the transport of chemotherapy to develop a comprehensive model based on the miscible displacement equations in the liver. We then apply the simulation to data from MD Anderson patients diagnosed with hepatocellular carcinoma who have undergone transarterial chemoembolization treatment. This final model enables us to provide insights into the evolving dynamics of TACE within the liver.