In recent earthquakes, tunnels in earthquake-prone zones have suffered significant damage. The cracks in these tunnels contain a mix of transverse, longitudinal and inclined cracks, suggesting a complex 3D tunnel response in these zones. In this article, we semi-analytically solve for the complete 3D response of a shallow lined tunnel embedded in a linear elastic homogeneous half-space under seismic waves that are incident on the tunnel from arbitrary directions. The waves scattered by the tunnel are represented using cylindrical waves, and the coefficients of these waves are estimated by enforcing the stress-free condition at the tunnel and ground surface, and the compatibility condition at the soil-liner interface. Accurate enforcement of the ground surface stress-free condition is accomplished by solving an improper contour integral that requires special treatment at branch points, removable singularities and poles. The convergence of the solution series is investigated, and the converged stresses and displacements at various locations in the soil and tunnel liner are compared with the results from past semi-analytical and numerical studies for 2D and 3D scenarios. Robustness and accuracy of the algorithm over a wide range of tunnel liner-to-soil stiffness ratios, tunnel depth-to-radius ratios, tunnel liner thickness-to-radius ratios, incident wave frequencies and directions are then tested to demonstrate the applicability of this algorithm as a fast and reliable tool for the initial design of a tunnel in an earthquake-prone zone.