In assessing risks of exposure to toxic contaminants in environmental systems, the transport dynamics of contaminants through complex heterogeneous media is a key determinant. Quantitative predictions of contaminant concentrations and fluxes are often required for systems involving different classes of physical dynamics and sharply discontinuous material properties. For many volatile contaminants, special conservation laws apply at interfaces between neighbouring system subdomains. These conservation laws couple dependent variables between the subdomains in a dynamic sense. Within the subdomains, essentially arbitrary physical processes may occur, possibly coupling many dependent variables. We describe the development of a code applicable to such multiphysics problems. The code uses an advanced black box solver as an engine to solve the partial differential equations appropriate to each subdomain. Simple geometric domain decomposition is used to define subdomains. In this first implementation, parallel instances of the solver engine (that is, in each subdomain) communicate synchronously with each other to facilitate time-stepping of the total system solution. The interfacial conservation laws act as constraints to the subdomain solutions, effectively providing dynamic internal boundary conditions. Convergence characteristics of the discretised interface algorithms are discussed in terms of validations against recently derived analytical solutions for diffusion-limited transport in partitioning laminates. Finally, an application of the code to coupled benzene vapour transport in variably saturated porous and permeable media is outlined.