Functional ultrasound (fUS) has emerged as a powerful neuroimaging modality due to its high resolution in both space and time, low cost and potential portability. Nevertheless, fUS signals provide only indirect observations of neuronal activity through the neurovascular coupling, and hence require the blind separation of latent neuronal sources while also deconvolving their hemodynamic responses. In this work, we propose a data-driven convolutive block-term tensor decomposition-based model for multi-trial fUS measurements, where each source has a spatiotemporal representation comprising a low-rank spatial map and a piecewise-constant neuronal activation signal convolved with a trial- and source-dependent hemodynamic response function (HRF) with a physiologically plausible shape. We propose a constrained optimization framework for the model computation, which consists of alternating projected gradient descent iterations. Simulation results are reported that demonstrate accurate recovery of spatial maps and reliable estimation of activation temporal profiles across various noise levels, while confirming that HRF estimation remains the most challenging part of the problem.