We derived working equations for calculating field-dependent NMR chemical shifts in paramagnetic molecules in solution, leading to two different computational approaches. The "finite field approach" is valid for arbitrarily high field strengths and arbitrarily low temperatures, but has the disadvantage that it requires numerical integration to account for the orientational averaging in solution. Furthermore, we also developed a "2nd order approach" that is only valid up to 2nd order in the external magnetic flux density B 0 , i.e., in the weak field / high-temperature regime. This approach has the advantage that the orientational averaging can be performed analytically. Furthermore, at this level of theory, the field-dependent part of the shift separates cleanly into two additive contributions: An "indirect effect" that is due to the partial orientation of molecules with an anisotropic susceptibility tensor in solution and which takes exactly the form that is well-known in the literature, and a "direct effect" that originates from the nonlinear response of the system to the external field and depends on a generalization of the fourth-order τ tensor introduced by Ramsey, Vaara, and others for diamagnetic molecules. In order to calculate this τ tensor, we derived analytical sum-over-states equations for fourth derivatives of the electronic Helmholtz free energy, which generalize the Van Vleck and Van den Heuvel-Soncini equations that are examples for second derivatives of the Helmholtz free energy. We expect the new sum-over-states equations to have applications also outside the realm of NMR. We tested the newly developed approaches using the NiSAL-HDPT complex as an example, describing the paramagnetic shifts as the sum of contact shifts (modeled within the spin Hamiltonian approximation) and pseudocontact shifts (modeled within the ligand field theory approximation and employing the point-dipole approximation). It turned out that the finite field approach converges very quickly with the size of the integration grid and already a modest 26-point Lebedev grid is sufficient for accurate results. Furthermore, we could show that at room temperature and at the field strengths that can routinely be achieved, the 2nd order approach is sufficient for describing the field dependence of the shifts. Finally, we could demonstrate good agreement between the experimental and calculated field dependence of the shifts in the NiSAL-HDPT complex.