Identification of a transfer function
The objective is to produce noise, presumably Gaussian distributed, that has both white and 1/f components. This is evident from the desired power transfer function, i.e., the square of the OP's "S" expression, which would be used to convert input white noise to output white plus 1/f noise. The  OP does not state whether the PSD should be unilateral or bilateral, so it will be assumed, without loss of generality, to be unilateral. Hence the unilateral noise PSD is

$$\mathrm{S^2 = k_1^2(1+k_2/f)} \tag{1}$$

where $\mathrm{k_1^2}$ is the unilateral noise PSD of white noise only, with units of $\mathrm{V^2/Hz}$ and $\mathrm{k_2}$ is the noise corner frequency, with units of $\mathrm{Hz}$. In a digital simulation,

$$\mathrm{k_1^2} = 2\sigma^2 /f_s = \sigma^2 /f_{nyquist} = 2\sigma^2 \times \Delta t    \tag{2}$$

where $\sigma$ is the standard deviation of the Gaussian white noise (units: V), $\mathrm {f_s}$ is the constant sampling frequency (units: Hz), $\mathrm {f_{nyquist}}$ is the Nyquist frequency (units: Hz), and $\mathrm {\Delta t}$ is the constant sample time spacing (units: s).

To generate the noise, first select parameters in accordance with equation 2. For example, if it is desired to have $\mathrm{k_1^2} = 5.00 \times 10^{-8}$ $\mathrm{V^2/Hz}$ and $\mathrm {f_s} = 1$ $\mathrm{kHz}$, then $\sigma = 0.005$ $\mathrm{V}$. This is the standard deviation entered into the Gaussian white noise generator in the figure below:



The second step is to specify $\mathrm{k_2}$. Using the noise generator described previously here, it is not feasible to directly enter a noise corner frequency into the non-white noise generator code. However, it is necessary to enter the desired low limit and high limit for the non-white noise, as shown in the dialog box in the figure, so the high limit can be adjusted empirically to obtain the desired noise corner frequency. Note that $\gamma = 1$ generates $\mathrm{1/f}$ noise and there is an optional gain adjustment as well.

The next figure shows a temporal trace of the noise produced with $\mathrm{k_1^2} = 5.00 \times 10^{-8}$ $\mathrm{V^2/Hz}$, $\mathrm {f_s} = 1$ $\mathrm{kHz}$, $\sigma = 0.005$ $\mathrm{V}$ and the non-white noise parameters shown in the first figure:



Then the third figure shows the average of 1000 32k FFT PSDs:



The noise corner frequency is where $\mathrm{f = k_2}$. From the plot, this is about 1.6 Hz or so.

If this noise was the dominant one in an experiment, and if the signal was baseband, then the measurement would be afflicted with $\mathrm{1/f}$ noise and all the annoying consequences that entails. But if it was feasible to modulate production of the signal, before the signal encountered the noise, then the signal could be deliberately located in the white noise 'flats' and a lock-in amplifier could perform the customary demodulation.