Background: Thermoembolization presents a unique treatment alternative for patients diagnosed with hepatocellular carcinoma. The approach delivers a reagent that undergoes an exothermic chemical reaction and combines the benefits of embolic as well as thermal-and chemical-ablative therapy modalities. The target tissue and vascular bed are subjected to simultaneous hyperthermia, ischemia, and chemical denaturation in a single procedure. To guide optimal delivery, we developed a mathematical model for understanding the competing diffusive and convective effects observed in thermoembolization delivery protocols. Methods: A mixture theory formulation was used to mathematically model thermoembolization as chemically reacting transport of an electrophile, dichloroacetyl chloride (DCACl), within porous living tissue. Mass and energy transport of each relevant constituent are considered. Specifically, DCACl is injected into the vessels and exothermically reacts with water in the blood or tissue to form dichloroacetic acid and hydrochloric acid. Neutralization reactions are assumed instantaneous in this approach. We validated the mathematical model predictions of temperature using MR thermometry of the thermoembolization procedure performed in ex vivo kidney. Results: Mathematical modeling predictions of tissue death were highly dependent on the vascular geometry, injection pressure, and intrinsic amount of exothermic energy released from the chemical species, and were able to recapitulate the temperature distributions observed in MR thermometry. Conclusion: These efforts present a first step toward formalizing a mathematical model for thermoembolization and are promising for providing insight for delivery protocol optimization. While our approach captured the observed experimental temperature measurements, larger-scale experimental validation is needed to prioritize additional model complexity and fidelity.