Reacting porous media are widely distributed in nature, but retain a similar underlying dynamics: an interconnected matrix changes its morphology due to the chemical reactions occurring within itself. This leads to both changes in the internal porosity and in characteristic dimensions. The purpose of this work is to encompass all these phenomena in a single computational framework, with a particular focus on biomass pyrolysis and oxidation. Current numerical approaches for biomass particle pyrolysis focus either on the solid particle evolution or on the surrounding gas-phase dynamics, neglecting the coupled interactions between the two. The first category is mainly composed of simple one-dimensional models, which rely on empirical correlation to estimate the interface exchange terms. The second category solves only the surrounding flow field, often relying on a fixed interface assumption. This work solves the major issue related to both, encompassing them in a single numerical model. This is achieved through the use of a Volume-of-Fluid approach, tracking the moving interface between the porous region and the surrounding environment using a single grid. The anisotropic nature of biomass is also accounted for in the resulting multidimensional framework. Detailed chemistry and reaction kinetics are necessary for accurate results, and their incorporation requires ad-hoc numerical strategies. A novel approach is introduced to account for both internal porosity changes and interface regression simultaneously. The proposed model is validated against experimental data gathered from centimetre-scale particles in terms of temperature profiles, mass and volume variations. Results showed good agreement between simulated and empirical results, considering also the high uncertainties and variabilities related to biomass modelling. The proposed formulation is independent of particle shape and number, and can help estimate sub-grid models for higher-order models, often used to simulate reactor-scale phenomena.