<p>Topology optimization for mass and heat transfer has been implemented in three dimensional domains or on two dimensional planes with the lack of extension to 2-manifolds representing the curved surfaces locally similar to two-dimensional Euclidean spaces. In order to enlarge the design space and increase the design freedom, this paper develops topology optimization on variable 2-manifolds for mass and heat transfer in volume flow, where the volume flow is the fluid flow in a three dimensional domain. In the developed topology optimization method, thin-wall patterns are defined on variable curved surfaces represented as implicit 2-manifolds within the three-dimensional domain, where the thin-wall patterns are the structural patterns in the mid-planes of wall-shaped structures with ignorable thickness. The implicit 2-manifolds are homeomorphously defined on preset base manifolds. Fiber bundles is used to describe a thin-wall pattern together with the implicit 2-manifold as an ensemble defined on the base manifold. The topology optimization method on variable 2-manifolds is developed to optimize the fiber bundles for mass and heat transfer in volume flow. It is implemented by using a mixed interfacial condition that combines no-jump and no-slip types. The mixed form is achieved by the interpolation between these two types of interfacial conditions, where the interpolation depends on the material density representing the thin-wall patterns. Two design variables are defined for the thin-wall patterns and the implicit 2-manifolds, respectively. They are regularized by two surface-PDE filters. Variation of the implicit 2-manifolds is controlled by introducing the variable magnitude to the surface-PDE filter. The topology optimization problems are analyzed by using the continuous adjoint method to derive the gradient information of the design objectives and constraints. They are then solved by using the gradient based iterative procedures numerically implemented based on the finite element method. In order to use linear finite elements and reduce the computational cost, the variational formulations of the governing equations are stabilized by using the Brezzi-Pitkäranta, Petrov-Galerkin and general least squares techniques. These methods are applied in the three-dimensional domains, which are deformed according to the implicit 2-manifolds and described by Laplace’s equation. The adjoint equations are derived for the stabilized variational formulations of the governing equations. In the numerical results, the effect of variable amplitude of the implicit 2-manifolds and that of the Reynolds number, Péclet number and pressure drop are investigated to demonstrate the increased design freedom and extended design space.</p>