This paper is devoted to the study of a sub-diffusion equation involving a Clarke subdifferential boundary condition. It describes transport of particles governed by the anomalous diffusion in media with boundary semipermeability. The weak formulation of the model problem results in a time fractional parabolic hemivariational inequality. We first construct an abstract hemivariational evolutionary inclusion and prove its unique solvability using a time-discretization approximation, known as the Rothe method. In addition, a numerical approach based on a finite difference scheme in time and finite dimensional approximation in space is proposed and analyzed for the abstract problem. These results are then applied to establish the convergence of the numerical solution of the model problem. Under appropriate regularity assumptions, an optimal order error estimate for the linear finite element method is derived. Some numerical examples are provided to support the theoretical results.