The short answer is yes, the HLL scheme and its variants to eliminate oscillations and/or enforce the proper entropy solution can be used to solve systems with source terms on the right hand side. However, the devil is in the details and there are several places where the solutions can go off the rails.
So the long answer is that how you discretize the source term in time will impact the quality of the solution. If the source terms are stiff, in other words if the ratio of the largest to smallest magnitude eigenvalues is large, then the limiting factor in choosing the time step may not be the eigenvalues of the flux jacobian like it is for the problem without the source term.
Stiffness is not the only factor though. If the real part of the eigenvalues of the source term jacobian are positive, then the source term is a source. If the real part is negative, then the source term is a sink. If you have multiple equations, it is possible that you have both sources and sinks depending on the equation.
For sink terms, explicit temporal integration is usually a bad idea because it requires small time steps for stability. Too large of a time step and it can over-shoot the zero point and give you negative quantities, which is generally an unrealizable solution. Therefore, implicit integration is preferred, which will always give a realizable solution and permits larger time steps (although if stiff, smaller time steps will be desired for accuracy).
For source terms, the opposite is true. Implicit time integration is not desirable, while explicit time integration is unconditionally stable. But you may want to take smaller time steps still to preserve accuracy.
For mixed problems, you can try to split the terms into sources and sinks and do a mixed integration (IMEX temporal scheme). Or just bite the bullet and run implicit for everything, or explicit with small time steps.
This is all independent of whether the spatial terms are explicit or implicit in time. In fact, it is independent of how you discretize the spatial terms entirely -- this is why HLL should work fine, along with any other spatial discretization scheme.
And you can choose to use a fractional step, or operator splitting, method to treat the temporal integrations differently. This will introduce a steady-state splitting error proportional to the eigenvalues of the product of the operators, unless the spatial and source operators commute (not likely for real problems). The splitting errors can be corrected if needed, but it takes more work.
All of these issues are rather complex and subtle. Many papers and theses have been written about how to do it correctly, including from me. What is best for your problem may not be what is best for another problem. Hopefully this has enough answer to give you a jumping off point for further reading about methods to solve your particular problem. If you need help with the analysis or algorithmic aspects, the questions will be better suited for SciComp.