To give a big picture view of how the GPS solution is determined, consider the following equation:
$\rho_i = \sqrt{(x_i-x_u)^2+(y_i-y_u)^2+(z_i-z_u)^2} +c\Delta t$
where $\rho$ is essentially a range/distance from the user to the GPS satellite, $x,y,z$ are position coordinates, the subscript $i$ indicates the particular satellite, $c$ is the speed of light, and $\Delta t$ is a time delay.
Assuming that you have knowledge of the GPS space vehicle (SV), the $x_i,y_i,z_i$ values are known from the satellite ephemeris (this can be obtained from publicly available data, and more accurate ephemeris can be obtained via more secure methods). There are now 4 unknowns, implying that we need 4 GPS SVs to solve for the user location $(x_u,y_u,z_u)$ and time delay. More SVs can be observed, and an over-determined solution can be found from various numerical methods (e.g., a least squares solution), or a best-4 SV solution can be employed.
The time delay is essentially in the $\Delta t$ term. Various errors can be accounted for by augmenting the system of equations to include, but in no way limited to, ionospheric & tropospheric delays, relativity effects, and clock errors present in the receiver.
A multitude of simple and complex differential methods exist to essentially exploit similar delays between measurements and remove them without even solving for them (e.g., differential GPS and real time kinematics
Here is a short paper that discusses the observation equations and, more specifically, the GPS signal and code-generation.